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1 . Introduction 

The Taylor series method [1] has long been regarded as 
an efficient procedure for solving systems of ordinary 
differential equations. Frequently, it is necessary to 
algebraically manipulate the differential system into an 
equivalent system. The Taylor coefficients for this modi- 
fied system may be simply written. However, the required 
modification is a tedious and error prone task for all but 
the simplest systems. For this reason, the Taylor series 
method has often been excluded by numerical analysts from 
consideration as a general purpose integrator. 

In Moore [2] , the procedures for recasting a system, 
which is reducible to a rational form, have been described 
in detail. Barton, Willers, and Zahar [3] describe techniques 
for automatic step length prediction, local error estimation, 
and for choosing the proper number of terms in the series. 

The authors also include a comparison of the Taylor series 
method with the fourth order Runge-Kutta method [4] and the 
Bulirsch-Stoer rational extrapolation method [5] . For a 
wide range of accuracy., it was found that the Bulirsch-Stoer 
method generally required five times the amount of computing 
time, and the factor for the Runge-Kutta method varied from 
five to one hundred. 

A method for the automatic reduction of arbitrary 
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differential systems is described in Barton, Willers, and 
Zahar [6] . Also presented is a procedure to generate the 
computer routine which evaluates the Taylor series coeffi* 
cients of the reduced system. The system reduction and pro- 
gram generation are analogous to the output from a compiler, 
and the differential equations and initial values are treated 
as simple language statements that are input to the "compiler". 
The particular implementation in [6] is an interactive pro- 
gram written for the Atlas 2 computer in Cambridge, England. 

The target language is Atlas machine language code. 

In this paper, an implementation allowing wider usage 
is presented. The "compiler" is written in PL/I, and the 
target language is Fortran IV. In Section II, the reduction 
of a differential system to rational form is described 
along with the procedures required for automatic numerical 
integration. In Section III, the Taylor series method is 
compared with the Bulirsch-Stoer method and with the 
Nordsieck version of the Adams predictor-corrector method 
[7] for a number of differential equations. 

In Section IV, algorithms using the Taylor series method 
to find the zeroes of a given differential equation and to 
evaluate partial derivatives are presented. 

Section V discusses the PL/I implementation of the 
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Barton et al . algorithm. Appendix A contains an annotated 
listing of the PL/I program which performs the reduction and 
code generation. Included in Appendix B are listings of the 
Fortran routines used by the Taylor series method. Finally, 
Appendix C has a compilation of all the recurrence formulas 
used to generate the Taylor coefficients for non-rational 
functions which may appear in the defining system of 
equations. 
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II. The Taylor Series Method 

Consider the following differential system 

^ = f(t,y), y (a) =a, a<t<b 

where the are rational functions. To apply the Taylor 
series method to this system, the Taylor coefficients for 
the expansion about the point t Q =a are computed. The de- 
pendent variables y^ are then evaluated at t=t-^, with 

y .(t 1 )=f (VV 3 


( 2 . 1 ) 


3=0 T 
dt D 


ji 


The value t is now replaced by t^ and the p: )cess repeated 
until the y^ at the value t-^=b are evaluated 

Initially, it may appear that the applicability of the 
method only to differential systems involving rational 
functions is a severe limitation on the usefulness of the 

mA't't'ir'tH . y _ fnnpf i c on r»1-i o -i -n /->/-> o ev v ,,r ' ^4-^, 

“■*“ —■ — — — n V — — f — »"*■ W k-; VAWX& VA U -I.U / ^ V kJ f ^ ^ V- W*V * f 'w 

solutions of rational differential systems. Consequently, 
a large class of solutions of non-rational differential 
systems have equivalent representations as solutions of 
rational differential systems. 

To illustrate this point, the function y satisfying the 
differential equation 

_ (a sin(y) +e cos (y) 


( 2 . 2 ) 


dt 


, y ( o) = o, o<t<7T 

2 


( 2 . 3 ) 
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may be written as the function in the system 


du. 


_L 

dt 

= u 4 +u 5 

du 2 

= u 3 (u 4 +u 5 ) 

dt 


du 0 


dt 

= -u 2 (u 4 +u 5 ) 

du. 


dtT 

= u 4 u 3 (u 4 +u 5 ) 

du. 


dt 

= _U 5 U 2 


(2.4) 


u (0) = [0 ,0 ,1 ,1 ,e] 


0 < t < ti_ 
2 

u, u. 


k 2 3 

where u 2 -sin(u^), u^-cosfu^), u^=e ,u^=e 

To obtain the canonical system equivalent to (2«4) , auxil- 
iary variables are introduced so that each equation in the 
canonical system represents a single operation of either 
addition, subtraction, multiplication, or division. Once the 
canonical system has been generated and the order of evalua- 
tion determined, it is a simple task for the computer to pro- 
duce the formulas for the coefficients. 

In the implementation of the method it is important to 
determine how to best evaluate expressions of the form 
’max 
3=0 


max (i) i 

w- 2 Yi 3 <V ( W 


(2.5) 


where y. ^ (t )= — 
o' . , 
J • 


& <v 

dt D 
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should be a 


Also, it is necessary to decide whether j max 

constant value over the interval of integration or whether 

j should be changed at each integration step. Other 
max 

questions involve the procedure for varying step length and 
the method of estimating local truncation error. 


It was found for a number of test differential equations, 
including those in Section III, that Horner’s method [8] for 
evaluating (2.5) proved to be the most efficient. Horner's 
method applied to (2.5) is given by 

a) y i (t 1 )=y i (: ’max ) (t Q ) 

b) y i (t 1 )=y i (j) (t 0 ) + ( V t o )y.(t 1 ) 

for ^max -1 ""' 0 

Relative error in the Taylor series solution is con- 
trolled by methods analogous to those commonly used for other 
discrete integrators. The interval length is varied from step 
to step in order to yield a local relative truncation error 
less than some preset error bound. The error term resulting 
from the truncation to J max terms of the Taylor series for 


( 2 . 


y. (t,) expanded about t is 
11 o 

y i ( ^max +1) (5) (t x -t ) ^max +1 


t <£<t, 
o— — 1 


Thus, a local relative error bound of E requires that the 
step length h=t^-t satisfy 


h w +1 ± E min 

i 


N. 

i 

y . ^ max + ^ ( t T 
2 1 o 


N. = 
r 


f W for y i (t o } “ 0 
1 


for y i (t Q ) = 0 


(2 



wh 

y 


e i + v^ries over the set of indices for which 


max 


. ( 3 + 1 ) 

(t Q )+0. If y^ (t^)=0 for all i, then h is 


i o' ■ - 1 • o 

set to step to the end of the range. 


For the differential equations considered in Section III, 
the fixed j which proved to be most efficient was equal to 

lUaX 

the number of significant decimal digits carried by the 
computer. This was also found to be true for the equations 
tested in [6] . For many problems where large functional 
changes occur over the integration interval , and computation 
time is critical, a variable j max may produce a very efficient 
procedure. For a further discussion of numerical integration 
methods which are optimized by changing the order at each 
step, see [9] and [10]. 
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III. Comparison 

An age old problem confronting numerical analysts is 
the generation of effective procedures for the comparison of 
computational methods. It is virtually impossible to include 
such characteristics as simplicity of method, implementation 
effort, reliability, and efficiency in a conclusive evaluation. 
Almost all comparisons of numerical integrators are made 
solely on the basis of efficiency - usually measured by the 
number of integration steps or the computer time required 
to obtain solutions of equal accuracy. 

With third generation machines, the concurrent execution 
of programs, and optimizing compilers, the computer time 
required for solution is subject to wide fluctuations. These 
fluctuations are often of the same order of magnitude as the 
computation times being measured. Also, during the computa- 
tion, there is an overhead charge incurred when index regis- 
ters are saved, arguments are passed, and loops are generated. 

Many implementations of a numerical algorithm will reduce 
the overhead at the expense of generality. It is unfair to 
compare on the basis of computer time, routines which differ 
in their implementation philosophy, because for the moderate 
sized problems generally used as test cases the overhead is 
often a significant portion of the computation time. Conse- 
quently, a less general, low overhead method may perform 
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competitively with a less efficiently programmed and more 
general method. 

In comparing the Taylor series method with other methods, 
significant factors such as the extra storage needed, the 
difficulty in learning to use the "compiler", and the effort 
in debugging the Taylor coefficient routine if the "compiler" 
malfunctions, are difficult to include. Further difficulties 
result because the Taylor series method integrates a differ- 
ent system of equations than do the usual methods. 

To eliminate implementation dependence from the estimate 
of a method's efficiency, each test problem was integrated 
to determine the number of derivative evaluations required 
for solution. This number should be approximately constant 
for a given method regardless of implementation. For each 
derivative evaluation routine, the number of machine (360/91) 
cycles required for one pass through this routine was deter- 
mined. Table III-A shows the number of cycles required for 
some typical operations. 


9 



TABLE III- A 


OPERATION 

NUMBER ( 

D.P. Load and Store 

0 

D.P. Add and Subtract 

2 

D.P. Multiply 

3 

D.P. Divide 

12 

F.P. Add, Subtract, Load 
and Store 

1 

F.P. Multiply 

11 

F.P. Divide 

37 

D.P. Sin and Cos 

217 

D.P. Exponential 

217 

D.P. Square Root 

133 

D.P. Power 

400 


F.P. = Integer Arithmetic 

D.P. = Double Precision Floating Point Arithmetic 


*Not including overlap or simultaneous operations. 
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The number of cycles required to pass arguments from the 
calling routine was not counted and the overlap or simul- 
taneous execution of operations was not considered. 

The methods selected for comparison were the Bulirsch- 
Stoer rational extrapolation and the Nordsieck version of 
the Adams predictor-corrector. Both of these methods re- 
quire a number of functional evaluations to obtain a starting 
step size which satisfies the accuracy condition. If the 
initial estimate for the step size is far off, the number of 
evaluations used in starting could be quite large. For the 
problems considered here, the number of evaluations required 
to start were not counted. The test problems are five repre- 
sentative non-trivial differential equations encountered in 
a computation laboratory: 

Problem 1. Bessel Function 

y" = Y(2/t 2 - 1) 

Y { 0) = 0 
Y* (0) = 0 
y" ( 0) = 2/3 
0 < t < 25 tt/4 

Solution: Y(t) = sin(t)/t - cos(t) = tj Q (t) 


11 



Problem 2 


Coulomb Function [11] 


Solution: 


nuuitjm o 


Solution : 


Problem 4 


Solution : 


Y (-1 + l/t)Y 

Y (0) =0 

v’ (0) = (n/ie^-i) ) 1/2 

o < t < 20 

Y (t) = F Q (l/2,t) 


Restricted 3 -body problem [5] 

x" - X + 2Y* - a'(X + g) - a(X - g’) 

Y = Y - 2X* - a*Y - aY 

X (0 ) = 1.2 

x' (0) = 0 

Y ( 0 ) = 0 

Y (0) = -1.04935750983 
g = 1/82.45, g ' = 1 - g 

a = g/((X - g ' ) 2 + Y 2 ) 3/2 , a' = g’/((X + g) 2 + Y 2 ) 3/2 
0 1 t 1 6.192169331396 

The given range for t is one period. 


Y = -Y + (1 + t) costte^) 

Y ( 0) = 0 
0 < t < 5 

Y ( t) = e“ t sin(te t ) 
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Problem 5. A stiff equation [12] 


X = -2000X + 1000Y + 1000 
Y = X - Y 
X(0) = 0 

Y ( 0) = 0 

1 £ t <_ 4 

Solution: X(t)=l+A^e X 1 ^tA^o 

Y(t)=l+B 1 e” Xlt +B 2 e’ A2t 
X i=+2000 . 5001 ... 

X 2 =+. 49987500 ... 

A-^=“ .49975000 ... 

A 2 =-. 50024999 ... 

B^t. 00024993746 ... 

B 2 =-l. 0002499 ... 


Table III-B lists some of the results of testing four 
of the five problems. The column labeled "error" refers to 
the relative error of the computed solution at the end of 
the interval. "Cycles" is the number of machine cycles re- 
quired for each evaluation. "DE" refers to the number of 
evaluations required to integrate over the given interval. 
Finally, the column labeled "R" contains the ratio 

-Comparison- Cvcles -Comparison. 

DE ' Method ' cycles ( Method ) 


“ < T * y ;r a ir 1 “> - <*=i~ < r * y ;s^r io, > 
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TABLE III-B 



TAYLOR SERIES 

METHOD 


PROBLEM 

ERROR 

CYCLES 

DE 

1 

1.2xl0" 10 

789 

15 

2 

4.1xl0~ 9 

736 

10 

3 

2.7xl0“ 9 

23769 

103 

4 

4 . 7xl0~ 9 

2524 

557 



NORDSIECK 

METHOD 



PROBLEM 

ERROR 

CYCLES 

DE 

R 

1 

7.9xl0" 10 

22 

661 

1.2 

2 

5.5xl0“ 9 

19 

741 

1.9 

3 

l.OxlO” 9 

349 

2340 

0.1 

4 

2*4x10 9 

445 

18374 

• 

GO 


BULIRSCH 

-STOER RATIONAL 

EXTRAPOLATION 

METHOD 


PROBLEM 

ERROR 

CYCLES 

DE 

R 

1 

1.7xl0“ 10 

22 

1288 

2.4 

2 

1 . 6xl0 _9 

19 

790 

2.0 

3 

1.0xl0~ 9 

349 

5769 

0.8 

4 

1.6xl0“ 9 

445 

6612 

2.1 
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where the comparison method is either Nordsieck or Bulirsch- 
Stoer . 

The Taylor series method is superior to Nordsieck and 
Bulirsch-Stoer on Problems 1, 2, and 4 and inferior on 
Problem 3. Other results, not presented here, show that the 
Nordsieck and Bulirsch-Stoer methods are very inefficient for 
Problem 5, while the Taylor series method handles this prob- 
lem well. 

Once the user masters the fairly simple art of setting 
up the input for program generation, he has an easy means for 
applying the Taylor series method. If greater efficiency is 
required, the program may be optimized by the user who has 
some knowledge of Fortran. On the other hand, if an error 
occurs, the program may be difficult to debug. Finally, 
it should be noted that there exists an important class of 
problems where no Taylor series method program can at 
present be generated. In general, however, the method is 
a valuable tool for solving many problems and is certainly 
worth trying. 
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IV. Finding Zeros; Partial Differentiation 


In this section two algorithms are presented which may 
be used in conjunction with the Taylor series method. The 
first algorithm finds the zeros of a function and the second 
algorithm is used to find the partial derivatives of a 
function of several variables. 

The method used to solve for the zeros of a function 
is that of series inversion. The relevant theorem is quoted 
here without proof [13] . 

Given the power series 

f=f +E a. (t-t ) k (4 ‘ 1) 

° k=l k ° 

with positive radius of convergence and a^=fO, then there 
exists a unique power series 

t=t o + L b k (f -v k (4 - 2 > 


with positive radius of convergence and such that the two 

series are inverses in sufficiently small neighborhoods of 

t and f and b. = 1/a. . 
o o 11 


To develop a recursion formula for the coefficients b^ 
in (4.2), solve for (f-f Q ) in (4.1) and substitute into (4.2), 
resulting in 


t-t =1 b, 
° k=i k 


3=1 J 


(4.3) 
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Letting 


2 c jk (t-t o ) ] = 
]=k 


Z a. (t-t ) 

j=l 3 ° 


(4.4) 


for k>l 

and interchanging the order of summation in (4.3) leads to 


00 • i 

t-t =2 (t-t) 3 t c it b k 

o j=1 o k=1 3 K K. 


(4.5) 


Equating powers of (t-t ) in (4.5) yields 
b l" 1/c ll 

^ -1 (4.6) 

for j>_2 „ 

Rewriting (4.4) in terms of previously computed coefficients, 
we find 


b.= (E c..b )/c.. 
3 k-'L D K k J J 


l c., (t-t ) “*= Z c . t , (t-t )^ Z a. (t-t ) 


j=k 3k 


=k-i j ' k_1 


(4.7) 


j=l 


00 00 

=E Z c 
r=l s=k-l 


r+s 


s,k-l r ' o 


a_ ( t-t_) 

for k>2 

Substituting j=r+s and interchanging the order of summation 
yields 


00 . 00 . j-k+l 

Z c..< t— t ) 3 =I (t-t ) 3 Z c. 
• 7 lk O . , o - 

j=k j-k r=l 


j-r,k-l r 
for k>2 


(4.8) 


Finally, equating powers of (t-t ) yields 
j -k+1 

c jk =2 =1 c j-r ,k-l a r for k L 2 

j>k 


(4.9) 
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Also, note that 


c j 1 *aj j>l . 


(4 


The following summarizes the algorithm to find t' such 


that f(t’)=0 when the a are known, t Q is given sufficiently 
close to t' , and f 0 -f(t Q ). 


1) 

C ll : 

2) 

C j 1 

3) 

C jk 

4) 

b .= 

J 


=a- 


b l ly/c ll 


: . , t a for 2<k< j 

3 -r ,k-l r J 


5) 


r=l 

Ll C 3k b k /c j j 

Repeat 2) thru 4) for j=2,3... 

oo 

b ( -f 

C r o k=l D k l r o' 


i* 


To illustrate the application of this method, the 
differential equation for the ninth degree Legendre poly- 
nomial was integrated and the zeros of the function computed 
by series inversion. The results were accurate to the re- 
quested precision. 


For the computation of the partial derivative of a 
function of several variables f(y-^, y 2 ,...,y n ) with respect 
to y. , the Taylor series coefficients for the differential 
system 



dt 



j=l, . . . ,n 


. 10 ) 


. 11 ) 


(4.12) 
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are evaluated along with the coefficients for the function 
f(t). The derivative of f with respect to t may be written 


df 

dt 


n 

Z 

s — 1 


3f 

3y 
J s 



(4.13) 


Substituting (4.12) into ( 4 . 13 ) , it is clear that the desired 
partial derivative is the first Taylor coefficient of f. 

This procedure may be applied to any number of functions and 
was used to evaluate the Jacobian of the system given in 
Problem 3. The results of this computation were as accurate 
as the input data. 
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V. PL/I Implementation 


The program to generate a Fortran subroutine which 
evaluates recursive Taylor series coefficients for a system 
of differential equations has been written in PL/I. The 
PL/I language was chosen, instead of a string processing 
language like SNOBOL, because PL/I contains an adequate 
set of string manipulating functions and because of the 
similarity between PL/I and Fortran statements. Since the 
PL/I statements are Fortran-like, changes may be incorpor- 
ated into the processing program to suit individual needs, 
with greater facility than might otherwise be the case. 

In the implementation of [3], the defining system may 
contain derivatives of arbitrary order and the differentia- 
tion operator may appear on the right hand side of the equa- 
tions. Without a serious loss in generality, the current 
implementation is restricted to systems of first order 
differential equations and the differentiation operator 
may not appear on the right hand side of the equations. 

The program reads in the defining system of equations 
from the PL/I SYSIN data set, and the equations are checked 
for balance with respect to parentheses, but no determina- 
tion is yet made as to whether they represent valid expressions. 
The program then attempts to generate the Fortran subroutine 
to evaluate the Taylor series coef f icients . 
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It is not a difficult task to add information about the 


interval of integration, the accuracy required per integration 
step, etc. to the input definition of the system of equations. 
This information may then be edited into appropriate driver 
routines. The next job steps are compilation and execution 
of the Fortran program. However, these additions are de- 
pendent upon the computer installation and upon individual 
requirements. The PL/I program is written and annotated so 
that modifications, similar to the ones just mentioned, are 
fairly straightforward. 

The first input card to the PL/I program is a control 
card containing the words DIFFERENTIAL EQUATIONS, which may 
appear anywhere in columns 1 to 72. Sequence numbers are 
permitted in columns 73-80. The word EQUATIONS may be 
optionally followed by the letters SP or DP, which is a 
request to generate a single or double precision routine. 

The default value is single precision. The nomenclature 
for the i— equation in the differential system is 
Y (1 , I) =f (T, Y) , where the first subscript of Y denotes 
differentiation with respect to the independent variable T. 
f(T,Y) represents a valid Fortran expression. If it is 
more convenient to specify the system with a different in- 
dependent and dependent variable, say R and V, then it is 
necessary to include (R,V) on the first control card. The 
differential equations are specified next with a free form 
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format in columns 1-72. Each differential equation is ended 
with a semi-colon, except for the last equation, which is 
terminated with a colon. Any equations which are useful in 
defining the differential system may be included and are 
order independent. The next card following the differential 
equations is a control card containing the words INITIAL 
VALUES. The initial values are then specified in the same 
manner as the differential equations. 

To illustrate a sample input to the processing program, 

consider Problem 3 above written as four first order equations. 

The input data required to specify the construction of a 

double precision routine may have the following form 

DIFFERENTIAL EQUATIONS DP(T,U) 

U(l,l) = U (3) ; U (1 , 2) = U (4 ) ; 

U ( 1 , 3 ) = U(l) + 2 . D0*U (4) - AP* (U (1) + G) - A*(U(1) - GP) ; 

U (1 , 4) = U (2) - 2 . D0*U (3) - U (2) * (AP + A); 

G = 1 .DO/82. 45D0; GP = 1.D0 - G; 

A = G/DSQRT ( (U (1) - GP)**2 + U(2)**2)**3; 

AP = GP/DSQRT ( (U ( 1 ) + G)**2 + U(2)**2)**3: 

INITIAL VALUES 

U(l) - 1.2D0; U ( 2 ) = 0 . DO ;U ( 3 ) = O.DO;U(4) = -1.04935750983 

The generated Fortran routine will have the structure 
SUBROUTINE COEFF ( U, ITSMAX) 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION U (ITSMAX, 1) 

Fortran statements necessary to compute the 
1 to ITSMAX Taylor coefficients for the equiva- 
lent canonical system given U(I) , 1=1,4. 

RETURN 
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ENTRY INITAL ( T , U , ITSMAX) 

• Initialization of all Taylor coefficients 
: to zero followed by the assignment of the 
| ^initial values specified as input data. 

RETURN 

END 

With the above routine and the two Fortran routines listed 
in Appendix B, edited in the appropriately indicated places, 
a complete Fortran program for the numerical integration of 
the sample problem may be developed. 

Standard output from the processing program contains 
listings of the defining differential equations and the gen- 
erated Fortran routine. Since it is necessary to have some 
measure of the computer time required per pass through the 
COEFF routine in order to properly assess the effectiveness 
of the Taylor series method compared to other popular methods 
for solving differential equations, an operations count in 
terms of additions and multiplications is also printed. 

The process by which the Fortran routine is generated 
is very similar to the way a compiler generates assembler 
language routines. For a complete description of the 
algorithm see [3] . The differential system is reduced to 
canonical form, which is the representation of the system 
in terms of the elementary operations of + - * /. The 
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decomposition is accomplished by the method of bounded con- 
text translation [14] . The next step consists of an elimin- 
ation of redundant operations from the canonical system. 
After the system has been optimized, a tree search is per- 
formed to determine the computational order. For some equa- 
tions, it may be desirable to examine a number of the inter- 
mediate quantities in this process. Coding DEBUG in the 
FARM field of the processing programs EXEC statement will 

produce this listing. For the Riccati equation. 

2 2 

y 1 = y + 3t , the DEBUG listing has the form given in 
Appendix D. 

RMAT is the procedure which performs the decomposition 
of the differential system. LEVEL denotes the current level 
of recursive calls to the procedure. The integer K denotes 
the element in the equation being scanned. TYPE is an 

j. t 

integer representing the K— element. Table V-A is a 

listing of the correspondence between the integers and the 

elements. The E in (C,0 ,v|e) denotes the print mode that 

lists the input equation to the procedure. The equation is 

enclosed by the delimiters #$. C,0,V represents the print 

mode that lists the K — element which is either a constant, 

th 

operator, or variable. If the K— element is a constant, 
it is replaced by # % . The integer i designates the position 
of this constant in a tabulation of all constants that 
appear in the differential system. The constant table is 
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TABLE V-A 


TYPE 

-1 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 


ELEMENT 

constant 

variable 

+ 


* 

/ 

( 

) 

# (left delimiter) 

$ (right delimiter) 

% (function specifi 
cation) 

* * 
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listed after the optimized canonical form. 


The heading on the right indicates entries into the 
recurrence matrix, where the canonical system is eventually 
stored. R represents the row of the matrix, OP the opera- 
tion, and A ( 1) , A(2) the two possible arguments. A(3) is the 
name associated with this operation. If there is no external 


name associated with this operation, the name is generally 
represented as ?r, where r indica Les the row in the matrix 
storing the result of the operation. The symbol $ in the 
recurrence matrix is used for the composite operation = „ 

The first differential equation processed is the one for the 
independent variable, which makes the system autonomous. 
After the last equation has been processed, the complete 


recurrence matrix is listed both before and after it has 


been optimized. As mentioned earlier, the constant table 
is listed at this point. 


The next step involves searching the recurrence matrix 
to initialize the matrix D described in [3] . The D matrix 
is used to determine the computational order of evaluation 
of the coefficients. The dimension of the matrix is the 
number of rows in the recurrence matrix. Briefly, starting 
with the result of the operation for a given row in the 
recurrence matrix, the arguments of the operation are 
traced backwards thru the recurrence matrix to ascertain 
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TABLE V-B 


% 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 


Function 

exp 

Iog 10 

log e 

sin 

cos 

tan 

sinh 

cosh 

tanh 

sqrt 
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their dependence upon other operations. The DMAT ENTRY 
statement lists the row currently being initialized, and 
finally the entire D matrix is listed. To aid in an inter- 
pretation of the recurrence matrix, a table is constructed 
showing the correspondence between this matrix and the set 
of dependent variables in the canonical system. An integer 
pair, i j , in this table, indicates that the result of the 
operation in the i— row of the recurrence matrix is the j — 
dependent variable. 

In the reduction to canonical form, special functions 
which appear may cause their defining differential equations 
to be appended to the differential system. In this imple- 
mentation, the special functions are left in the reduced 
system and the corresponding coefficients for these 
functions are hard coded in the program generating routine. 

The symbol %j, where j represents an integer constant, is 
used to represent functions in the recurrence matrix. 

Table V-B shows the correspondence between the integers j 
and the functions they represent. 

This completes the description of the intermediate 
quantities required in the Fortran COEFF routine construction. 
The listings should be useful in debugging any malfunctioning 
of the processing program for a given differential system. 
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APPENDIX A 


taylcr: ppoc(parm) options* main); 

/******************************************************************* ** 


* ORIVER PROCEDURE USED TO CONSTRUCT A FORTRAN SUBROUTINE WHICH * 

* EVALUATES THE RECURSIVE TAYLOR COEFFICIENTS DERIVED FROM A * 

* SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS • FOR A DESCRIPTION OF * 

* THE ALGORITHM SEE • THE AUTOMATIC SOLUTION OF SYSTEMS OF * 

* ORDINARY DIFFERENTIAL EQUATIONS BY THE TAYLOR SERIES METHOD*. * 

* BY D. BARTON ET AL* COMPUTER JOURNAL V 14 (1971 > PP . 243-248 * 


OOOCOlOO 

0000020© 

©0000300 

00000400 

00000500 

00000600 

00000700 

00000800 


********* ************************************** ** ******************** /©©0OQ9PO 


OCL 

BFDC ENTRY (BIN FIXED) RETURNS ( CHAR( 15 1 VARI. 

BREA KF E NTR Y ( CHAR < *) VAR • CHAR ( * ) VAR* CHAR ( 4 1 VARl. 

CODE ENTRY(CHAR< *> V AR .BIT ( 1 ) , 8 IT ( 1 > ) . 

count entry < char < 4) var.chari*) var) returnsibin fixed). 

SPAN ENTRY! CHAR ( 4 ) var.char<*» VAR»CHAR<4) VAR.CHAR(*> VAR): 
DCL 

(CS.WS) CHAR (400 ) VAR EXT ,OVRBL CHAR( 4 ) VAR EXT, 

IVR3L CHAR (4 ) VAR EXT.R( 500,4) CHAR(IS) VAR, 

RMAX BIN FIXED I NIT (O ) , 0< * .* ) BIN FIXED CTL. KFMAX EXT INITIO), 
ERROR BIN FIXED, CC BIN FIXED INtTUl.IED EXT, 

DEBUG BlT(i) EXT.NEQ EXT lNlT(0),KO BIT! 1 ) .NSGMA EXT INIT(O). 
FL FILE OUTPUT.SN CHAR(4) VAR.CB CHAR ( 15 ) VAR; 

DCL 

(NMUL.NADD, NETS, NMTL, NATS, NATL) INIT(O) EXT. PARM CHAR(IOO) VAR 
LBLM3) LABEL I NI T (LW ,LW. LD) • LBLB (3 ) LABEL I NIT( L X.LD.LO ) * 
CST(IOO) CHAR! 25 ) VAP EXT.IC EXT IN IT < 1 ) . NEQTNS EXT; 

/* 

CALL stint; 

IF INDEX! PARM, 'DEBUG* ) -i=0 THEN OEBUG=*I*B; 

/+ READ IN THE SYSTEM OF EQUATIONS 4/ 

csTd > = *o,5* ; 

CALL I NPUT! ERROR, CC ) « 

IF IED=0 THEN CST<1 )=»0.5* ; ELSE CST! 1 )=• O. 500* ; 

GO TO LBL A(ERPOR) ; 

LWt NEQTNS=COUNT(CS.* Y<1. • )*i; 

/4 INSERT DIFFERENTIAL EQUATION FOR THE INDEPENDENT VARIABLE TO MAKE 
THE SYSTEM AUTONOMOUS 4/ 

CS=»Y(1 , * | | BFDC (NEQTNS > | I • 1=1 #0 • M * $ • | | CS | | • A • 1 
/4 READ IN THE INITIAL VALUES 4/ 

cc=2; 

CALL INPUT(ERROR,CCI ; 

CS=CS | |DVRBL| 1 * (* ft BFDC (NEQTNS ) ( f • )=• f f IVRBLl | • :* ; 

GO TO LBLB! ERROR) ; 

lx: ipass=i ; 

/4 FACTOR EACH DIFFERENTIAL EQUATION INTO ELEMENTARY OPERATIONS 4/ 
LY: DO WHILF(SUBSTR!CS.1,1 )-*=•#•»? 

)€Q=NEQ+t ; 

CALL 0REAKF(CS. • , WS); 

CALL RMAT! •#* f fwsf f • 4* . R.» MAX, KO ) l 


OOOCIOOO 
00001100 
00001200 
00001300 
00001400 
00001500 
00001600 
00001700 
C0001800 
00001 90 0 
00002000 
00002100 
0000220© 
00002300 
. 00002400 

00002500 
00002600 
4/00002700 
00002600 
00002900 
00003000 
00003100 
00003200 
00CC330C 
00003400 
00003500 
00003600 
00003700 
00003600 
00003900 
00004900 
00004100 
00004200 
00004300 
00004400 
00004500 
00004600 
00004700 
000C4600 
00004900 
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IF NEQ=1 THEN R<RMAX.4)=IVRBLI 
END ; 

IF -.DEBUG THEN GO TO LZ5 

PUT EDITC •RECURRENCE MATRIX* I < SK IP< l ) * X( 60 ), A > J 
DO 1=1 TO RMA X * 

PUT EDI TCI »R(I*l)*R<l*2)«R<I«3t«R(l*4) ) <SKIP,X<60), 

F<2> * X<1> *A<2) • X<5) »3 A(l5)i; 

END ; 

/* eliminate redundant operations from the recurrence matrix */ 

LZ5 CALL OPTMZE <R *RMAX) ; 

if -tDEbug then go to te; 

PUT EDIT! ‘OPTIMIZED RECURRENCE MATRIX*) C SK IP ( 1 I • Xt 60 > ♦ A ) I 

la: no i=i to rmax; 

PUT ED IT < I *R< I ♦1)*R(I*2)*R<I*3)*R(I«4)) (SKIP* X( 60) *F( 2) • 

XC 1 ) ,A<2> *X<5 ) *3 A 05) ) » 

end : 

lb: if -.debug then go to lbc; 

PUT EDITt 'CONSTANT TABLE*) <SKIP(2>«A); 

DO 1=1 TO ic; 

PUT EDIT< * A * *1 »* = * * CST <1)1 <COLUMN<MOO< I-t *4) *290 ) * 

A .F(2 ) .A. A) : 

END ; 

LBC : ALLOCATE O < RM AX *R MAX ) : 

/* GENERATE THE MATRIX O WHICH IS USED TO DETERMINE THE ORDER IN WHICH 
THE TAYLOR COEFFICIENTS ARE COMPUTED */ 

CALL OMAT t R * R MAX * D ) * 

IF -.DEBUG THEN GO TO lc; 

PUT EDIT< *0 MATRIX*) <SKIP<2)*A); 

PUT FOIT< <<I , ’, ** J*0( I , J) 00 TO RM AX) 00 1=1 TO RMAX) ) 

(SKI P,8 <F< 2 )«A.F<2) *F<41 ,X<5) ))$ 

/* GENERATE THE FORTRAN ROUTINE TO COMPUTE THE TAYLOR COEFFICIENTS 
CALL C OGE < R *R MAX »D * KO) ? 

IF -.KO THEN GO TO LO? 

LC : FREE D ? 

IP*SS*:IPASSTl * 

IF IPASS>2 THEN GO TO LOT 
CALL CCDE < *C* • *0*B* *1 * B) I 

CALL CODE < * ENTRY I NITAL < • II IVR9L t I • * • 1 1 DVRBL |J * * ITSMA X) * . 

*0*B*«1*B> 5 

CALL CCDE < * 00 2001 ITS=1 • ITS MAX* * *0* B** l*B)» 

CALL CODE < * 00 2001 I XV*I ♦ * f | BFDC< RM AX )• • 0* B* • 1 « B )? 

CALL CODE < *2001 Y< ITS *TXV=0.0 • **0 *8* * l *B )? 

DO MH I LE < C S-»= * * ) i 

CALL BPEAKF < CS * • * * * MS ) * 

CALL SPANf fcS* *('•*)* «SN) * 

CALL 8REAKF < MS * * 1 * * CB) * 

CALL CCDE (DVRBL | t # < 1 * • I | SN M • ) ■ H WS* * I • B* ' 0* 8 > 5 
END 5 

CALL CODE (DVRBL | f • ( 2 * • | f EFDC< NEQTNS ) | | • )=1 *0* **0 , B**0 , B); 
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CALL 'CODE < * RETURN* ,»©*B, • I*B|J 

CALL CODEC* END* , *0 • 9* * I • B K 

PUT EDIT( * OPE RATI ON COUNT <OC> FOR ONE PASS THRU THE CuEFF * 

H 'ROUTINE*) <SKIP(3)«A); 

PUT EDITC *ITS*AX - THE NUMBER OF COEFFICIENTS COMPUTED*) 

<SKIP<2> »A); 

PUT SDIT( *AS - AN AOOITION OR SUBTRACT ION • > <SKIP*A)‘. 

PUT EDIT! *MD - A MULTIPLICATION OR DIVISION*) <SKlP*A>; 

WS^*OC - < • MBFOCINADD) It * ♦ < • I 1 BFDCf NATL ) | | • 4 - * \ | BFOC( NATS ) f | 
**ITSMAX)*ITSMAX/2 l*AS ♦ < • II BFOC CNMUL ) I I • 4- < • |) BFDC C NMTL ) I I 
* 4 - • | |BFDCC NMTS) ff * *TTSH AX > * ITS MAX/24MD* ; 

PUTEDITOlS) ( SKI P ( 2 ) * A) * 

ld: end taylor; 


00069900 

00010000 

00010100 

00010200 

00010300 

00010400 

©OOIOSOO 

00010600 

00010700 

00010690 

00010900 

00011000 

©9011100 
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* process; 

TAYLOR: PROCIPARM) OPT IONS ! MAIN) ? ©OOOOIOO 

/********* 444 * 4 * 44 *****44*4 4*4* 4444444*4*******44* **4*4 44* *4444**44* *4 ©0000200 

* DRIVER PROCEDURE USED TO CONSTRUCT A FORTRAN SUBROUTINE WHICH * ©0000300 

* EVALUATES THE RECURSIVE TAYLOR COEFFICIENTS DERIVED FROM A * 0000040© 

* SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS* FOR A DESCRIPTION OF 4 000005*© 

* THE ALGORITHM SEE * THE AUTOMATIC SOLUTION OF SYSTEMS OF * 0CCOO60C 

* ORDINARY DIFFERENTIAL EQUATIONS BY THE TAYLOR SERIES METHOD • * * ©000070© 

* BY D. BARTON ET AL. CCMPUTER JOURNAL V 14 11971) PP. 243-248 4 0000080C 

*4 *4 44 44* 4* ***** * ********* 4* *4 444 4* *4* 4*4 * 4*4 444* ******************** /O00OO9OO 


C CL ir W ^ v * w % 

B*OZ ENTRY! BIN FIXED) RETURNS < CHAR! 15 ) VAR)* ©C0C110© 

BREAKF ENTRY(CHAR< 4) VAR*CHAR<4) VAR. CHARI*) VAR I* 00001200 

CODE ENTRY! CHARI *1 VAR *BIT! 1 ) * BIT ! 1 ) ) • 00001300 

COUNT ENTRY! CHAR !*) VAR*CHAR<4) VAR) RETURNS I BIN FIXED)* 00001400 

SPAN ENTRY! CHARI 4) VAR.CHAR!*) VAR*CHAR!4) VAR«CHAR!4) VAR); OOOCISOO 

DCL 00001600 

!:S*WS) CHAR (409 I VAR EXT *DVRBL CHAR<4) VAR EXT* ©0001700 

IVRBL CH AR < 4 ) VAR EXT*R< 509*4) CHAR! 1 51 VAR* 00001800 

RMAX BIN FIXED I N1T!0 ) .0! *.* ) BIN FIXED CTL* KFMAX EXT INITIO)* 00001900 
ERROR BIN FI XED*CC BIN FIXED INIT!I)*IED EXT* 00002000 

DEBUG BIT!!) EXT *NEQ EXT INITCO)*KO BIT(1)«NSGMA EXT INITIO)* 00002100 
FL p ILE OUTP UT *SN CHAR!4) VAR*CB CHAR! 15) VAR « 00002200 

DCL 00002300 

!NMUL*NADD.NMTS*NMTL* NATS* NATL) INITIO) EXT* PARM CHARI100) VAR* 00002400 
LBLAI3I LABEL INIT ILW *LW*LD) *LBLBI3) LABEL I NITI L X* LD ,LD ) . 00002500 

CSTII00I CHAR! 25 ) VAR EXT*IC EXT IN IT ( 1 * , NEQTNS EXT; 00002600 

/* 4/00002700 

CALL STINT; 00002800 

IF INDEX (PARM* 'DEBUG •)-.«© THEN 0E8UG=M*8; 00002900 


/* READ IN THE SYSTEM OF EQUATIONS 4/ 
cstii )= , e*5* ; 

CALL INPUTIERRCR«CC) ; 

IF IED=0 THEN CST (l ) = *0*5 # * ELSE CST 1 1 )=• 0 • 5D0* 5 
GO TO LBLAIERROR) ; 

LWt NEQTNS^COUNT ( CS * , Y!1** ) *1 ; 

/* INSERT DIFFERENTIAL EQUATION FOR THE INDEPENDENT VARIABLE TO MAKE 
THE SYSTEM AUTONOMOUS */ 

CS= »Y!I ** | fBFDCINEQTNS ) [ | • )-l *0 * | f * ? • | |CS| I'M*? 

/* READ IN THE INITIAL VALUES */ 

cc=2; 

CALL I NPUT (ERROR* CC ) * 

CS=CS | |DVRBL l| * < * If BFOCINEQTNS H | • )*• | \ IVRBLl | * 5 • ! 

GO TO LBLB(ERROR) ; 
lx: i pa 55=1 ; 

/* FACTOR EACH DIFFERENTIAL EQUATION INTO ELEMENTARY OPERATIONS 4/ 
LYt DO WHILE! SUBSTR(CS*1*1 A* I ; 

NEQ = NEC*1 ; 

CALL BREAKF (CS ** ;*.WS>; 


WVVJUMV 

00003100 

00003200 

00003300 

00003400 

00003500 

00003600 

©0003700 

00003800 

©0003900 

©0004000 

00004100 

00004200 

00004300 

©0004400 

00004500 

00004600 

00004700 

00004800 
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CALL RMATC I !*Sl 1 • * * « R*RMAX* KOI ; 
IF NE 0=1 THEN R< RMAX ,4 ) = 1 VRBL? 


end; 

IF -.DEBUG THEN GO TO LZ; 

PUT EDIT! 'RECURRENCE MATRIX*) < SK IP< 1 ) , X< 60 ) * A > ? 

DO 1=1 TO RMAX; 

PUT EDITH ,R<!.n.R(I.2),RU.3),R(!.4)) CSKIP*X<60>« 
F(21,Xlt) *AC2) *XC5>*3 AC15)); 


END ; 

/* ELIMINATE REDUNDANT OPERATIONS FROM THE RECURRENCE MATRIX */ 
LZ5 CALL OPTM 2E C R * R MA XI* 

IF -*DEBUG THEN GO TO LB; 

PUT EDIT! 'OPTIMIZED RECURRENCE MATRIX* I (SKIP! 1)» X(60)«A) I 


LAI DO 1=1 TO RMAX? 

PUT EDITH *R(I ,1 >*RC1«2)*R<1 *3)«RCI*4> ) C$K IP* Xf 60 > *F< 2) . 
XC1) *AC2) *XI5> ,3 A<15)H 

End ; 

lb: if -oebug then go to Lee; 

PUT EDIT! 'CONSTANT TABLE*) <SKIPC2)*A>5 

do 1 = 1 to ic; 

PUT EDIT! •#•*! *•=• *CST CI>) I COLUMNCMODC I-I • 4)*29+l ) « 
A*F<2)*A*AI ; 


END ; 

lbc: allocate ocrmax*rmaxi ; 

/* GENERATE THE MATRIX D WHICH IS USED TO DETERMIME THE ORDER IN WHICH 
THE TAYLOR COEFFICIENTS ARE COMPUTED */ 

CALL OMATCR*RMAX*D) ; 

IF -.DEBUG THEN GO TO LC: 

PUT EDIT! *D MATRIX*) <SKIPC2)*A); 

PUT EDIT! CC1«****J«0CI*J) DO J=1 TO RMAX) DO 1=1 TO RMAX)) 

C SKIP* 8 <F<2)*A*F<2>*F<4I*X<5»>>; 

/* GENERATE THE FORTRAN ROUTINE TO COMPUTE THE TAYLOR COEFFICIENTS */ 
CALL CCG6 C R *R MAX *0 • KOY ; 

IF -.KO THEN GC TO LD! 

lc: free d; 

IPASS=1PASS+1 ; 

IF IPA SS>2 THEN GO TO LD; 

CALL CODE I *C* **0*0* *1 ■ B) ; 

CALL CODE < * ENTRY I NI TAL I * 11 I VRBL H • * • M DVRBL f f * * ITSMA XI • • 


•0*b • * i * bi ; 

CALL CODE C • 

CALL CODEC * 
CALL CODEC *2901 


OO 2001 ITS=1 * ITS MAX ***0*B** 1 • B ) • 

DO 2001 IXV=1 , • H BFDCCRMAX)* *0*B* *1*0); 
VC ITS * I XV=0 *0 • ♦*0*B. * 1 *b); 


DO WHILECCS-.= * * ); 

CALL 8REAKF CCS* • « * « MS ) * 

CALL SPAN! WS»* C * ••) • «SN ) * 

CALL SREAKFCWS**)* *C0i; 

CALL CODEC DVRBL 1 1 • C 1 * • \\ SN (1 * ) * I I WS. • 1 • B* * O* 0 ) ; 


END ; 


00004000 
00005000 
00005100 
00005200 
00005300 
00005400 
©0005500 
00005600 
00005700 
00005800 
00005900 
00006000 
00006100 
00006200 
00006300 
00006400 
00006500 
00006600 
00006700 
f 0906800 
00006900 
00007009 
00007100 
00007200 
00007300 
00007409 
00007500 
00907600 
OOC 0770© 

00007800 

00007900 

00008000 

00008109 

90008200 

00008300 

00008400 

00008500 

00008600 

00008709 

©0008899 

00098909 

00009000 

0000919© 

00009200 

00909309 

00909400 

00009509 

00009699 

©0009700 
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CALL CODECDVRBL! | M2* • | | BFDC ( NE CTNS ! | | • 1= I *0' . • 0* B . * 0 ■ 0 ) ; 

CALL CODEC • RETUR N * • * © * 8 * * 1 » 8 1 ; 

CALL CODEC* END* * *0 • B * * 1 * B > ; 

PUT EDITC * OPE RATI CN COUNT (OC> FOR ONE PASS THRU THE COEFF * 

| J *ROUTI NE * > CSKIP(3),A>: 

PUT EDITC *ITSNAX - THE NUMBER OF COEFFICIENTS COMPUTED*! 

< SKI P(2> ,A> ; 

PUT EDIT( *AS - AN ADDITION OR SUBTRACT ION * J (SKIP* A) I 

PUT EDIT< *MD - A MULTIPLICATION OR DIVISION*) (SKIP, A); 

*S=»OC = ( • || BFDC(NADD)| f * ♦ ( * I I 8FDC ( NATL ) I | • ♦ • | f BFDC < NATS I | | 
•*ITSMAX)*ITSMAX/2 )*AS «■ ( • | | BFDC ( NMUL ) II * ♦ < • 1 | BFDC ( NMTL > f | 

5 « j fSFDCI NMTSli | rSWAXj *1 TSMAX/2*MD * $ 

PUT EDIT(WS) (SKIP(2),A); 

ld: end taylor; 


oootseoe 

© 000990 © 

00010000 

0001010C 

0001020 © 

0001030 © 

© 001040 © 

00010500 

©0010600 

00010700 

00010800 

00010900 

00011000 

©0011100 
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* process; 

code; proc < string. sum. cwmnt) recursive; 

/* PROCEDURE TRANSFORMS THE INPUT * STRING* INTO FORTRAN CARO IMAGES */ 
DCL 

EXTRACT ENTRY! CHAR <*) VAR. CHAR I*) VAR* CHAR (A I VAR)* 

LlBF ENTRY(CHARI *) VAR. BIN FIXED). 

SIGMA ENTRY! CHAR <*) VAR). 

OC1_ ST CHAR (4 09 ) VAR. <SUM*CMMNT > BIT ( 1 ) .STR IMG CHARI* > VAR* 

SN CHARI 7) VAR.FL FILE OUTPUT EXT.IED EXT 5 
DCL NSGMA STATIC BIN FIXED EXT. SYSA 8 IT <1 I EXT.IC INITIl). 

SON BIN FIXED STATIC INITI19000). EISI3) CHARI8) VAR EXT; 

/* */ 
IF -.SYSA | LENGTH! STRING XII THEN GO TO LA? 

SN= SUB STR I STRING. 7. 5) ? 

IF SN=*SUBRO* ! SN= * BLOCK* THEN PUT PAGE; 

la; if -.sum then go to lc; 

lb; call extract i stri ng *ei s i 3 ) *st ) ; 

IF ST*** THEN GO TO LC 5 
CALL SIGMAIST); 

GO TO lb; 

LCI IF -£MMNT then go to ld; 

SON=SQN+100 ; 

PUT FILEIFL) EDITISTRING.'OOD^SQNICSKIP.A.COLUMNITSI.A.FIS)); 

IF SYSA THEN PUT ED IT C STRI NG. • 900 ■ .SON )<SK I P. A.COLUMNI 73 ) • A .FI S ) > 

return; 

ld ; sn= * • ; 

IF -CMMNT THEN CALL Li BE ISTRI NG* IEC) ; 

00 1 = 1 TO 6 WHILE CVERIFYISUBSTRCSTRING. I. I). *0123436799 • )=9) ;END 
IF I-.=l THEN SN=SUBSTR|SU8STR|STRING. 1. 1-1) II SN.1.6); 
string=substr<string.i > ; 

LS=L£NGTH I STRING) ; 

DO 1*1 TO LS BY 65; 

ST=SN || SUBSTRI STR ING. IC.MINC65* LS — IC*1))* 

SQN=SQNMOO; 

PUT FILEIFL) EDIT I ST . • 009 • • SON > C SK IP * A • COLUMN I 73 ) * A . FI 5 ) ) ; 

IF SYSA THEN PUT EDIT 1ST . • 900* . SON ) (SKI P. A .COLUMN I 73 ) . A. FI 5 ) ) 


00000190 

00000290 

00000300 

00000490 

00000 soo 
00000600 
00000790 
00000600 
00000990 
00001000 
00901190 
00001200 
00001300 
00 001400 
00001590 
00001600 
00001790 
00001690 
00001900 
00002000 
00002100 
00002200 
00002300 
100002400 
00002500 
00002600 
00002700 
;00 002809 
00002900 
00003000 
00003100 
00003200 
00003300 
00003409 
00003500 
;00CQ3690 


IC=IC«-65; 

IF 1=1 THEN SN = • X 

END ; 

END CODE ; 


00003700 

0000380Q 

00003900 

00004000 
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* process; oofooi^o 

coge: proc<r*rmax*d*koi ; ©00002©© 

/* PROCEDURE GENERATES THE FORTRAN TAYLOR COEFFICIENT ROUTINE */ 00600300 

OCL OO(M1Q40O 

SIP DC ENTRY! 8 IN FIXED) RETURNS ! CHAR! IS ) VAR I, ©0600560 

BFTC ENTRY! BIN FLOAT 1 53 ) • BIN FIXED! RETURNS ! CHAR ! 50 ) VAR1* 00000600 

BREA KF ENTRYICHAR!*) VAR * CHAR ! I ) • CHAR ! * I VAR** ©0000700 

CODE ENTRY!CHAR< *1 VAR *BIT< 1 > * BIT ! 1 ) I « ©06068O© 

COLA ENTRY RETURNS <BI N FIXED)* ©6CC096O 

FUDGE ENTRY RETURNS! BIT !1 M * ©000106© 

OP E NTRY < CHAR < 1 ) ) RETURNS (BIN FIXEO)* ©©©©*1©© 

RiPLACE ENTRYICHAR!*) VAR*CHAR !* ) VAR*CHAR<*> VAR*BIT!1))« 00001200 

SPAN ENTRYICHAR! *1 VAR. CHAR!*) VAR.CHARC*) VAR* CHAR!* ) VAR I * 00001300 

OCL 00001400 

R!*.*> CHAR!*) VAR *D! *.*) BIN FIXED. DR!*) BIT!1) CTL « KO BI T< 1 ) . 00601560 
CTBL ! *) CHAR !15) VAR CTL. DVR 9L CHAR ! 4 ) VAR EXT. I DMAX* RMA X) 00001600 

BIN FIXED*!*!*) *C(*) ) BIN FIXED CTL . FTBL ! 100 * 2 ) BIN FIXED EXT. 00001700 
L3LI11) LABEL I NIT 4PMS. PMS. MPY ,DVD. EOL.ERR . 00001800 

ERR. ERR . I NT . FN.E XP ) *CST ! 109 I CHAR 425 * VAR EXT. CO CHAR ! 1 ) * I ED EXT .00001000 
L(4) LABEL. DC! 19 ) CHAR41) EXT.WS CHAP! 400) VAR E XT* NEQTNS EXT* 00002000 
DCL 00002100 

LEPMAX<3> I NIT43.4.2) .LEP INIT!0).CH CHAR! II* 00002200 

IVFN CHAR! 10 ) VAR* I VRBL CHAR !4 I VAR EXT. 00002300 

DEBUG BI T< 1 ) EXT *OT ! *.* ) BIN FIXED CTL. 00002400 

CAB CHAR !1) EXT •! NNUL .NADD.NMT $♦ N MTL* NAT S. NATL ) EXT. ©0002500 

IQ BITfl) INIT! *1 *B) * ADS CHARISO) VAR.EPLG CHAR! 200) VAR? ©00 0 2600 

DCL 00002700 

( ARG ! 2 ) CHAR !25l VAR, ! CA3 * ZERO ) BIT! 1 ) * ( LHS. LHSARG ) CHAR! 15 ) VAR .00002800 
!LP< 2) *RP<2) *UO! 2 ) ) CHAR!!) V AR. K CMAX * KFMA X* T5S CHAR(5) VAR 6©9©29©0 

) FXT; ©0003000 

f 0003100 


ZERO=«l *8 ; 

EPLG= 

•SUBROUTINE COEFF! * || DVP8L II ' . ITS MAX ) ; IMPL IC IT REAL* • I f C48 1 | 

M A-H.O-Z) ;• 11 'DIMENSION • j j DVRBL f | • < ITS* AX, 1)1 1000 DO 2000 Ml 
•its=2*itsmax;itsmi=its-i ;itspi=itsm;fitsmi=float!Itsmi j; • 1 1 
•2000ccntinue;return; • ; 

1 VFN=DVRBL| | • !• I | BFDC! NEQTNS) | f ■ )• % 

dmax=dimi o.n ; 

ALLOCATE DR <0 MAX) *CTBL I RMAX ) , M ! DMAX ) * C! OMAX )* DT ( RMAX* RMA X) ; 

dr=*o«b; dt=o; 

/* CORRESPONDENCE T AQ_E BETWEEN R!I*4) & Y«J) */ 

KC = NE QTNS-1 ? 

DO 1 = 1 TO RMAX * 

IF SUB STR ! R f I * 4 ) *1 . MlNlLENGTHfR < I , 4) l*LENGTH!DVRBL ) ))=OVRBL 
THEN CALL SPAN ! R f I *4 )*•<•*•)•« CT BL ! I ) K 

else oo; kc^kcm; ctbl<ii=bf©c<kc>; end; 

end ; 

IF DEBUG 


©06032© 0 
©0003360 
06603460 
©0003500 
00003600 
00003700 
©0003800 
00003000 
©0004066 
00604100 
00004200 
00004360 
00004400 
©0004560 
06004660 
©0004760 
06004860 
©0004060 


A- 8 



THEN DOS PUT EDIT! 'CORRESPONDENCE BETWEEN RECURRENCE MATRIX ROWS 'COCOS©©© 
| | • A NO THEY ARRAY') <SKIPC2).A>; ©©CCS1©© 

PUT EDIT I ( I *CT0LI I I DO 1*1 TO RMAX)) <SKIP*12 < F< 3 > ♦ XC 1 ) *©06052©© 
A<3)*X<41)>; 000653©© 

End; ©©0654©© 

/* EVALUATION OF THE SET M */ ©066550© 

KM=0; ©066560© 

DO J=l TO DMAX; ©0©C57©6 

DO 1=1 TO OMAX; ©0665800 

IF DP ( I » THEN GO TO LB; ©CC6590© 

IF OT<I*JI>2 THEN GO TO LC; 00CC6000 

lb: end; ©cocci©© 

GO TO LE* ©000C20© 

LC: KM=KM+i; MfKM»=j; ©00C6306 

LE: END! ©0C66406 

IF KM=0 THEN GO TO LEA; O0OO65OC 

PUT EDIT! •** PROLOG IS NOT CURRENTLY IMPLEMENTED **' > I SK IP ( 2 ) * A) ;©0C C 6606 
KO='0*B; 0C0C67C© 

return; fococeo© 

lea: PUT PAGE EDIT < • ** LISTING OF THE GENERATED FDRTRAN ROUTINE **•> 00066900 

IA>: PUT SKIP; ©QC67000 

lf: lep=lep+i ; oocc7i©c 

DO 1=1 TO lepmaxilepi; 00CC72*6 

CALL BREAKFCEPLG.* ; **WS ) ; ©066730C 

CALL CCDEOlS*'© '8* '©'8) ; ©00C7400 

END ; ©0007500 

IF EPLG-*= f * THEN GO TO LG1 ©0007600 

FREE DR.CTBL*W*CiOT; 00067700 

RETURN; ©00C7800 

/* EVALUATION OF THE SET C */ 00CC7900 

lg: kc=o: ©ocosoo© 


IF ZERO THEN TSS='<t*'5 ELSE TSS=MITS**; 

DO 1=1 TO DMAX ; 

IF DR( 1 I THEN GO TO LI 1 
I s -ZERO 6 10 

then do; if r<i*i> = *s* then oo; KCMAX=i; GO to lk; end; 
GO TO Li; 

eno; 

DO J=l TO DMAX ; 

IF DT<I*J>>© THEN GO TO LIT 
end; 

kc=kc+i; kcmax=i ; CIKCI*!;- 
li : end : 

IF TO S— ZERO THEN DOT I0='0'B; GO TO LG* END* 

IF KO0 THEN GO TO LK; 

PUT EDIT! • ** EQUATIONS ARE NOT WELL POSED *♦« 1 (SKIP(2)*A); 

ko='©«b; 

PUT EDIT! # D MATRIX*! <SKIP€2)*A); 

PUT EDIT((U.SM.DUiJ) DO J=1 TO RMAX) DO 1*1 TO RMAX I ) 


coooeio© 

00668200 

00CC83©© 

09CC84C0 

OOCC8SO© 

©ooceco© 

ooooezo© 

006C88©© 

©0668900 

©0CC90C6 

00009100 

00009200 

COC09300 

00669406 

C00695QQ 

C0009CCC 

COC097C© 

OOOC980© 
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lk: 


ll: 


lu ) 

L<2> 
L< 3) 
L34S 

L<4) 

fn: 
ERR : 


< SKI P .8 CFC2 ) • A. Ft 21 .FC4) »X< 5) > l* 

IF DEBUG 

THEN PUT EDITCCCl *• .*. J.OTfl^JI DC J=1 TO RHAX) DO 1= 1 TO RMAX1 ) 
< SKIP. 8 CFC2).A.F<2) ,FC4).XC5)li; 

return; 

LHSARG=CTBLCC0L4<RCKCMAX.4I.R.RMAX >15 
LHS=DVRBL | |TSS| f LHS ARG | | • )=• ; 

CA3 = < SUBSTRCR CKCMAX .3 ) . 1 * 1 » = • A* >; 
ka*o: UO.LP.RP*** ; 

IF RCKCfAX. !»=•*** 

then do; nop-is: cas« •- end; 

ELSE oc: NCP=CP?R ?KC“AX.I5 ? ; Cu=OC{NuP»; End; 

DO 1-1 to 2 ; 

IF RCKCMAX.l >-*X* e 1-1 THEN GO TO LL; 

IF SUBSTR<R<KCMAX, IM 1. 1.1 ! = •#• 
then do; ka=ka+i; 

IF NOP <3 S -.ZERO 

then do; argci>-**; if nop*i|i *2 then co-**; end; 
ELSE arg <il*CST< substrcr ckcmax.h-u* 2) >; 

end; 

ELSE do; if rikcmax. IH1 )»IVFN 

THEN ARG < I > — CTBL ( C0L4 1 1 VRBL* R. RN AX > I; 

ELSF ARG < I >=CTBL<C0L4IR<KCMAX, 1+ 1 > . R , RMA X 1 1 ; 

ch^substr crckcmax. i-m i.i, i >; 

IF CHst*t | CH=* — * 

then do; uoci »sch; lpci»=*c*; rpcd=*j*; end; 

end; 

end ; 

IF KA=0 THEN KA-4; 

I F NOP=9 THEN GO TO INT; 

IF NOP— 1 1 THEN GO TO EXP; 

IF NO P-5 


000099PC 
DOC 10000 
© 001010 © 
00010200 
00010300 
00010400 
00010500 
00010600 
00010700 

oooioeoo 

©001090© 
©001100© 
00011100 
00011200 
00011300 
00011400 
00011500 
00011600 
00011700 
00011800 
00011900 
©0012000 
00012100 
0001220 C 
00012300 
©0012400 
00012500 
00012600 
00012700 
00012800 
00012900 
00013000 


THEN DC; IF -.ZERO THEN GO TO LT ; LH5=* • ; LP < 2 > .RP < 2>= * • ; END; ©001310© 

GO TO L<KA»; 00013200 

SIF NOP— 4 ©0013300 

THEN WS®*-SIGFA<I xv*2. ITS; • | IDVRBLII • C I XV * • t I ARG( 2 1 | | • 1** | | DVRBL 1 100013400 
•ITSP1-IXV, • 1 ILHSARG 1 !'))>/'! iDVRBLl | •( 1, • | | ARG< 2 ) 1 | • >•; ©0013500 

ELSE WS-ARGC 1 ) f |CO t |LP<2) f |UO<2) | |OVRBL| |TSS | ) ARGC 2) | |* J • | |RP< 2) ;©0C136ft© 
GO TO L34; ©0013700 

; WSs=UO<l > | IDVRBLl 1TSS | | ARG 1 1 ) | f * )* IfCOf |LP<2> M ARGC 2) | |«P<2) ; 00013800 

GO TO L341 00013900 

: IF -.ZERO THEN GO TO LT ; 00014000 

WS»LHSl |U0C1» I 1 APG(1> | | CO | |LP(2)| |UO<2)| I ARGC 2 1 1 | RP C 2 ) * 00014100 

IF NOP <3 THEN NADD»NArO+i; 00014200 

IF NOP<5 e NOP>2 THEN N*UL=NMUL+1 ; 00014300 

GO TO LS; 00014400 

: GO TO LBLCNOP1 ; 00014500 

IF -.F UOGE ( CTBL *DR .FTBL *R *RMAX » OT J THEN RETURN ; ELSE GO TO LU1 00014600 

PUT EDITC ILLEGAL OPERATOR IN COGE ) <SKIP<2)*A>; 00014700 
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KO**0*B; 00014800 

RETURN; 00014960 

✓» EQUALITY OF TWO SERIES */ 00015000 

EQLt WS=UOC2> | lOVRBLt 1TSS M ARG<21 f | * >• ; ©0015100 

GO TO LS? 00015200 

/* ADDITION OR SUBTRACTION OF TWO SERIES */ 00015300 

PMS: WS*U0C1» ! |DVR8L| |TSS| I ARG<1>| 1 • » • 1 f 00| |LP«2)| |UO<2>| IOVRSLI |TSS|| 00 01540© 
ARG< 2) | i • f |PP< 21 ; ©0015500 

NADD«NADD*1 ; 0001560© 

GO TO LS? 00015700 

/* MULTIPLICATION CF TWO SERIES */ 00015800 

MPY5 IF -.ZERO 0001590© 

then do; ws=* sigma; ixv*i .its; * I IdvrblI I mixv, 'IIargci HI ' )*• 1 1 oooieooc 

DVRBL || M ITSP1 —I XV * * j | ARGf 2 ) | |« I 1*5 00016100 

nmts=nmts«-i ; NWTL*NMTL«-1 ; 00016200 

END* 00016300 

ELSE 00016400 


MDZ: DO; WS-UC< I) | I DVRSL ( f TSS I f ARGO ) I f * I * ( (CO | | LP<2 J | f UO( 21 f | 

DVRBL | |TSS| | ARGf 2)| j • >• | |RP<2>; 

NMUL=NMUL+1 ; 

end; 

GO TO ls; 

/* DIVISION OF ONE SERIES BY ANOTHER */ 

DVDS IF ZERO THEN GO TO WDZ; 

WS=*< • | |DVRBL | | MITS* * |1 ARG<1 > | | 

•)-SIGMAMXV=2.ITS;« | lOVRBLl | » <IXV.* | | ARG(2l| |« l**| I 
3VRBL | | MITSPl-I XV* • | 1 lHSARG| 1 *>))/• | | DVRBLi j Ml* • j | ARG< 2) j | • 1 •; 
NMTS=NMTSO 5 NMTL-NMTL-1 7 
GO TO ls; 

I NT: IF ZERO THEN GO TO LT ; 

IF SUBSTR < R < KCMAX *31*1 * 1 1* * THEN GO TO LT* 

WS=U0<2> | lOVRBLl | M ITSMl *• f |ARG<2) | | * ) /FITS Ml • ; 

CALL CODE<LHS | |WS**1 # B**0*B>; 

go to lt; 

/+ SERIES RAX SEO TO A POWER */ 

EXPt IF KA“3 

THEN DO! IF -.ZERO THEN GO TO LT ; 

WS— LPO ) I |UO(l ) I I ARG < 1 ) | | Rpf 1)| I 1 |LP«2I| I AftG<2)| | 

rp<2> ; 

NADO=NADD«-20O ; 

GO TC ls; 

end; 

IF ZERO 

THEN DO: WS«LP<11 | |UO< 1 I 1 | DVRBL | 1TSS| | ARG ( 1 1 | 1 • ) • | | RP < II | | 

•*** | |LP(2» | | ARG <2 > 1 | RP(2I ; 

NADD — NADD4-200 ; 

end; 

ELSE DO? IF IED=1 THEN CALL REPLACE C ARG <2 ) . • D* • « E» • • 0« B I ; 

ADS^BFTC <ARG<2) fl *0. IEDM 


©0016500 
00016600 
©00 1670 Q 
00016800 
0001690© 
00017000 
00017100 
©0017200 
00© 17300 
©001740© 
0001750© 
00017600 
00017700 
00C17 80 0 
00017900 
OOO180OO 
00018100 
0001 820© 
0001830© 
0001840© 
00018500 
00018600 
00018700 
00018800 
0001890© 
00019000 
0001910© 
00019200 
0001930© 
000 19400 
00019500 
00019600 
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IF 1ED=1 THEN CALL REPLACE! ADS** E* • • 0* , • 0 • B I J 
IF VERIFY! SUBSTR! ARG <2 ) *1 ,1 >*«♦-• >=0 
THEN ADS = * I* | | AOS j | • ) • ; 

*S^* SIGMA! IX V=1 • I TSMl ; <» | ] ARC! 2 ) | | •-! IXV-1 )*» \ \ AOS f | 
•✓FITSMt ) *• | IDVPBLI | • (IXV,* | | LHS ARG | | • )*• | | DVRBL | | 
MITSP1 -I XV,* | | ARG! 1 ) | | • >1/* | lOVRBLl f Ml** f 1 ARG! I) | |« ) 
NATS=NATS+35 NATL— NATL— 3 J NMTS=NMTS-*-3? NMTL=NMTL-3J 

end; 

LS: CALL CODE (LHS || VS , *1* B,*0« Bi; 

lt: dpc kcmax* = *i *e; 

00 1 = 1 TO DMA X ’ 

IF -OR(II C OT(I ,KCMAX)>^D THEN OT ( I *KCM AX ) = DT( I ,KCMA X>— 1 ? 

end: 

LU: DO 1 = 1 TO OMAx; 

IF -OR! I) THEN GC TO LG? 

END ; 

zero=*o*b; 
dr=*o*b; dt=d; 
go to lf; 

END COGE ; 


00(11*700 

00019800 

00019900 

00020000 

00020100 

100020200 

00020300 

00020400 

00020500 

00020600 

00020700 

00020800 

00020900 

00021000 

00021100 

00021200 

00021300 

00021400 

00021500 

00021600 
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* process; 

C0L4: PROC<CS*R*.RMAX> RETURNS ( SIN FIXED!? 

/ 44 ******** ******************************************* ********* ****** 

* THE PROCEDURE SEARCHES COLUMN 4 OF THE RECURRENCE MATRIX TO FIND * 

* THE ROW NUMBER WHICH CONTAINS CS * 


©OO0O10C 

00000200 

00000300 

OO0O04*« 

00000500 


44444*44* 4444 44*44444 *4444 4**4 444444**444444*4444 4444444444444444444/ 0000060-0 


DCL CS CHAR(*» VA R • I . RM < 4, * > CHAR<*» V AR* RM AX BIN FIXED-. 

CT CHAR< 151 VAR? 

CT=SUBSTR ( CS « 2— VERI FV ( SUBSTR ( CS. 1 * t!« 1 | ; 

DO 1=1 TO RMAX? 

IF RM!I*4)=CT THEN RETURN! I )? 

END ; 

PUT EDIT! *44 ERROR IN INPUT EQUATIONS • * • | | CS || 1 ’ * CAN NOT BE 
•FOUND IN COLUMN 4 OF THE RECURRENCE MATRIX 44*) <SKIP(2!«A>? 
STOP? 

END C0L4? 


000007*0 
00000800 
OO0OO90O 
00001000 
00001100 
OOCC120C 
U ©0001300 
00001400 
00001 500 
00001600 
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* process; ooeoeioe 

omat: procirm,rmax,dm> ; 00000200 

/ ** 4 * ************ ******* ************ ********************************** C0C0030C 


* PROCEDURE CONSTRUCTS THE KATRIX D WHICH IS USED TO DETERMINE ♦ 00000400 

♦ THE ORDER OF TAYLCR COEFFICIENT EVALUATION * O0O005"© 

****************************** ************************************* **/OOQ0©6OC 

DCL ©0000700 

BFDC ENTRY! SIN FIXED) RETURNS I CHAR ( 1 5 I VAR I* ©OfOCeOO 

COLA ENTRY RETURNS < 01 N FIXED)* 00000900 


/* 


LA: 


l„C ? 


S a AN ENTRY(CHAR( ♦) VAR *CHAR ( * ) VAR»CHAR I * ) VAR,CHAR<*> VAR ) 5 
DCL 

RM(*.*1 CHAR! *) Y*R S EFAX SIN FIXED , t BIN FIXED, 

R BIN FI XED,L(50O,2> BIN FIXED, CH CHAR <1 ), II VRBL, DVR BL ) 

CHAR C 4 ) VAR EXT,JCH CFARI2) VAR,VN CHARI 4 ) VAR, DEBUG BITH) EXT 


©0O01000 
OOOflllOO 
©0CC120C 
00001300 
; 00001400 
4/00001500 


DW=-1 ; 

DO R=I TO PMAX; 

IF DEBUG 

THEN DO? IF R=l 

THEN PUT EDIT I ’ DMAT ENTRY *,R) ( SKIP! 2),A«F<2))« 
ELSE PUT EDIT 1 11 , * , R ) <A(1),F<3>); 

end; 
n=o ; i — l ; 

LI 1 , t ) -R * LU«2>=1 + (RMILU , 1 ), 1 )=« = » | RMI LI 1 ,1 )• 1 >=•%* 

I RMILU *1 I ,1) =• *' » ; 

IF RMILIt • 1) ,1 )-• » * • THEN N=N+1 ; 

L1=LII*1>; L2=L( I * 2 ) *1 ; 

CH=SUBSTR(RM(L1 ,L2)*1 *1 ); 

I* CH=«7» 

then do; lu*i ,h-substr(rmili,l2i,2>; 

IC4-COL4CPMIL1 ,L21 ,PM,RMAX); 

IF D“ {" * IC4 i\N+i l Htw DM t R » IC4 ) =N M ; 

LII+l ,21*1 -MRMILI I *1,1 ), 1 >=•=•! 

PMiLiiAi ,i >,i >=• x* | pm(h i+i, l ) ,i >=♦*» i; 
i *im ; 
go to la; 

end; 

CALL SPAN! M*| | RMCL1 ,L2 >, * ) • , • I * • VN I * 

IF VN-.s« * 

THEN IF VERIFYISUBSTRIVN.l ,1 )»0 THEN VN= SUBSTRI VN, 2 ) « 

IF VN=DVRBL 

THEN do; IC4=C 0L4IRMIL1 ,L2 ) ,RM,RMAX ) ; 

IF DM I R , I C4 ) <N THEN DMtR,IC4)«N; 

go to lb; 

end; 

IF CH= • A • THEN GC TO LB; 

IF RM(Ll «L2) *1 VREL 

THEN do; IF DMIR,1XN THEN om<r,1)=n; 

go to lb; 


©0001 600 
00001700 
00001600 
00001900 
C 0002000 
00002100 
00002200 
00002300 
00002400 
00002500 
00002600 
00002700 
00002800 
00002900 
00003000 
00003100 
00003200 
00003300 
00903400 
©0003500 
00003600 
00003700 
00003800 
00003900 
00004000 
00004100 
00004200 
©0004300 
©0004400 
00004500 
00004600 
00004700 
00004800 
00004900 
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fkd; 

ELSE DO; IC4=C0L4<RMfLl*L21*RM*RMAX>; 

L(IM *1 >-ica; 

GO TO lc; 

end; 

lb; if LU.2) = 2 ! R*«R. !>=•=• I FM< R, l ) = *5£ • 

THEN do; I =i-t ; 

I F 1=0 THEN GO TO LR; ELSE GO TO LBI 

end; 

L< I *2)=2; 

GO TO la; 

LR : END ; 

ND ; END D MAT ; 


00003000 

000051 oe 

00005200 

00005300 

00003400 

00005500 

00005600 

00005700 

OOOC50OO 

00005900 

00006000 

00006100 

00006200 
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* process; 

extract: proc t stri ng,wopd *exs ) ; 

/* PROCEDURE EXTRACTS THE SUMMATION OPERATOR PROM THE INPUT STRING 
DCL 

BF DC ENTRYtBIN PIKED) RETURNS I CHAR ( 15.1 VAR), 


BLNCO ENTRY C CHAR I* ) VAR) RETURNS! BIN FIXED), 

REPLACE ENTRY! CHARI*) VAR, CHARI*) VAR, CHAR I* I VAR*BITC1I>; 
DCL { STRING, WORD) CHARI*) VAR,EXS CHARI*) VAR, 

INSGMA,NSMTR ) STATIC BIN FIXED FXT,EIS!3) CHARI 8 ) VAR EXT; 


eooootor 

00000200 

♦/00000300 

00000400 

00000500 

00000600 

00000700 

©0000000 

00000900 


/ * 


T/ VVVVt WWW 


Exs=* • ; 

CALL EXT! WORD ,WRKRA,NRKRB) ; 

IF MR KRA=0 THEN RETURN; 

DO 1=1 TO 3; 

IF E IS 1 1 )= WORD THEN GO TO LA* 

CALL E XT I WORD , Ml ,LI > * 

ip mkmrkra then return; 
la : end ; 

E XS=S UB STR 1 STRI NG * MRKF A « MRKR8-MRKR A*1 ) « 

IF VOR0=EISI3) 

THEN CALC REPLACE (STRI NG.EXS, • SGMA • 1 | BFOC! NSGMA+1 ), «0 * B) ? 

IF WORO=EIS!2) 

THEN CALL REPLACE! STRI NG,EXS, • SMTP ■ | | 0FDCI NSMTRM )• *0 • B) 5 
IF WORD=E I St 1 ) THEN CALL REPLACEISTRlNGwEXSfl 0*B ) ; 

EXS=SUBSTRIEXS,INDEXIEXS*« < • ) *1 )5 
EXS=SUBSTRIEXS,1,LENGTHIEXS)-1 ) | | • ;* ; 

/* 

ext; proci w,mr ,id ; 

DCL V CHAR!*) VAR,MR,IL; 

hr s i l=© ; 

ld: if mrm <lengthistring) 

THEN I X=I NOEX < SUBSTR! STRING, MR* 1 )• V)« ELSE IX=0; 

IF IX=0 
THEN 

lda: do; mr=o; return; end; 

else mr=ix*mr; 

I L= MR — 1 ; 

le: ix=index! substristring,il*d, • >• ); 

IF IX=0 THEN GO TO LDA; ELSE IL=IL*IX; 

IF BLNCO! SU8STR!STRING,MR, IL-MR + 1 ) )-,=0 THEN GO TO LE; 

IF INDEX! SUBSTRtSTRlNG, MR, IL), •=« 1=0 THEN GO TO LD* 

IF MR =1 | VERIFY! SUBSTRtSTRlNG, MR-1, I ),•♦-*/)«=* ) = 0 THEN RETURN 

GO TO ld; 

END ext; 

END extract; 


oooouoo 
00001200 
00001300 
00001400 
00001500 
00001600 
00001700 
000C1800 
00001900 
©0002000 
00002100 
00002200 
00002300 
00002400 
00002500 
00002600 
*/00002700 
©0002800 
00002900 
00003000 
00003100 
00003200 
00003300 
00003400 
00003500 
00003600 
00003700 
00003800 
00003900 
00004000 
00004100 
; 00004200 

00004300 
00004400 
00004500 
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♦ process: 

FUDGE: PR0C!CTBL.DP.FTBL.R.RMAX.DT » RETURNS* BIT< 1 ) ) ; 

/** * ****** **** ***************** ************************************** 

* PROCEDURE CONSTRUCTS THE RECURSIVE TAYLOR COEFFICIENTS FOR * 

* FUNCTIONS SUCH AS EXP* SIN. COS. TAN ETC WHICH MAY APPEAR * 

* IN THE DIFFERENTIAL EQUATIONS * 

****************************** **************************************/ 

DCL 

CODE ENTRY! CHARI ♦ ) VAR • BIT !1 > . EIT<1 ) ) • 

COL4 ENTRY RETURNS I BIN FIXED): 

DCL 

DR ! * ) BIT(l) .CTBL!*) CHAR !* > VAR. RC* • * ) CHAR!*) VAR. FTBL < * * * » 

BIN FIXED. FN(I3.2> CHAP!6) VAR EXT.RMAX BIN FIXED. 

(ARG<2) CHAR !25) VAR. !CA3. ZERO) BIT<1). f LHS.LHSARG) CHAR! IS) 
V*R, ! LP ! 2 ) * RP! 2 ) . UO! 2 ) > CHAR! I) VAR. KCMAX. KFMAX. TSS CHAR ! 5 ) 
VAR. !I VRBL.OVRBL) CHAR14) V AR ) EXT.VNLHS CHAR C 15) VAR. 

LHSAI3) CHAR (15) VAR. FNL < 9) LABEL IN IT! EX* . L 19 . LN, 
SCT.SCT.SCT.HSCT.HSCT.HSCTI.PW CHAR!1) VAR* 

(NMUL.NADD.NWTS.NMTL. NATS. NATL) EXT.IED EXT* DT ! * . * ) B IN FIXED. 
ALPHA ! 2 ) CHAR129) VAR l NI T (• 4 . 34P945E- 1 *.* 4. 34294481 903251 80- 1 « ) 
/* *S 

DO KF=1 TO KFMAX; 

IF KCMAX<»FTBL!KF.2> 6 KCW AX>=FTBL!KF. 1 ) THEN GO TO LA • 

END ; 

PUT EDIT!*** FUNCTION NUMBER *. KCMAX. • IS NOT IN TABLE **•> 

! SKI P ! 2 ) *A*F (2) .A) ? 

RETURN! *0 *B) ; 
la : ; 

DO K“F TBL ! KF « 1 ) TO FT EL! KF . 2) * 

DR (K)= *1 *8 S 
DO J~l TO FM AX • 

IF -.DR! J) E. DT!J.K>>=0 THEN DT 1 J. K )=DT! J .K ) — l 5 

fnd; 

i=i+i ; 

LHSA!I )=CTBL!COL4!R!K.4 I.R.RMAX )>; 

END ; 

NF=R( KCMAX .2 ) l 
IF -.ZERO THEN GO TO LB 1 
1=0 ; 

DO K=FTSL!KF.n TO FTBLIKF .2)1 

1 = 1 + 1 : 

KX=R!K.2> ; 

VNLHS=DVR8l! ! * !1 .* I ILHSAII ) | | ■ )=* 1 
IF CA3 

THEN CALL CODE ! VNLHS f |FN!KX.21 1 !* (• I |UO!2) | | ARG! 2) | | ■ ) I ** *0*0. 

•o*8): 

ELSE CALL CODE ! VNLHS | |FN(KX.2) || • ! • ||UOf 2) 1 1 DVRBLM •! 1. • | | 

ARG ! 2 ) } |«l)*.*0*B.«O*B); 

END 5 


00000100 
88000200 
00000300 
00000400 
00000500 
00000600 
00000700 
00000800 
00000900 
00001 000 
0000110© 
000012*6 
00001300 
00001400 
00001500 
00001600 
0000170© 
0006180© 
00001900 
; 00©02000 
00662106 
00002206 
00002300 
00002400 
00002506 
©0002600 
00602700 
©0062800 
00002900 
00003000 
00063106 
©0003200 
00603300 
00003406 
00003506 
00003606 
00663700 
00063800 
00003900 
00004000 
000041*6 
00004206 
00604300 
00064406 
©0664506 
00664600 
00064706 
00004806 
00664906 


A- 17 



GO TO LZI 

LB5 IF CA3 THEN RETURN! *1*61? ELSE GO TO FNL(NF); 

/* SERIES COEFFICIENTS FOR EXPONENTIAL FUNCTION *f 
exp: call ccde<lhs 1 1 »sigma( ixv= 2 ,its;* I I* <ixv-t j*» I | 

dvrbl || mxv, *1 |arg<2> M * ) *• ! I dvrbl? | « ( IT SP1— I XV# • | | 

LHSA! 1 ) || M /FITSM1 * • *1*B,*9*B15 

NMUL=NMUL+1 ? KWTS=NMTS+2 5 NMT L=KMT L— 1 ; 

GO TO LZ? 

/* SERIES COEFFICIENTS FOR LCG BASE 10 FUNCTION */ 

LlOt CALL CODE ! • I p (ITS^EO.SI '?! LHS |? ALPHA! lED^l }'?** 5 \ ' DVR5L ‘ i ; ( 2 * 
I 1 AR G ( 2 ) M • ) / ' f | DVRBL 1 I • ! 1 • • | I ARG C2 )| | 

CALL CODE ( • IF tITS.GT.21 • f ?LHS 11 • <• ? j ALPHA! I ED4-I 1 1) * * • 
f | DVRBL | | TSS f ?ARG!2) ? | * I —SIGMA < IX V=2 , ITS Ml » I IXV-1 )** | ? 

DVR3L | ? MITSP1-I XV,* | j ARG <21 | ! *1** | | DVRBL f 1* Cl XV, • ? j 
LHSA! 1) II *1 ) /FITSMl > / • ? | OVRBL If * (1, * | | ARG!2 >| |* 1 • , • t *B, *0*8); 

NMUL=NMUL + 3; NAOD=NADD*l 5 NMTS=NMTSF2; NMTL— N MTL— 2 • 

GO TO lz; 

/* SERIES COEFFICIENTS FOR LCG BASE E FUNCTION *✓ 

LNt CALLCGDE!'IF (ITS. EQ. 21 • | ! LHS | I OVRBL ? I TSS \ \ ARG! 2 ) | I * >/ • 

| | DVRBL | | »<t , • | | ARG! 2 1 | | , i*»*0*P**9*Bl? 

CALL CODE < • IF (ITS.GT.2I • | |LHS | | • !• | | DVRBL I I TSS| | ARG! 21 

| | *1-* | | • SIGMA (I XV=2,ITSMi;( IXV-1 )** j | DVRBL? | •(ITSPI-IXV. * | ? 
ARG! 2 J | | •) *• | lOVRBLf | MIXV,* | | LHS ARG j j • ) 1/FITSM1/* | | DVRBL | j M 1 , 
| I ARG! 2) I f *1 /* I | DVRBL | I Ml,* ? ? ARG <21 j j • > * , • 1 • B, • 0 • B 1 ; 
NMUL=NMUL+2? NAOD^NADD 4-1 ; NMT SsNMTS+2 ? NMTL— N MTL — 2 » 

GO TO lz; 

/* SERIES COEFFICIENTS FOR THE SIN, COS, TAN FUNCTIONS */ 

sct: pm= «- * ; 

GO TO LC5 

/* SERIES COEFFICIENTS FOR THE HYPER SOL IC 3 INH, C3SH, TANH FUnCT IONS 

hsct: pm=**; 

LC: CALL C CK>E< OVRBL ! ?TSS |? LHSA! 1 1 ? | * 1-SIGMA! IXV=2, ITS? ! I XV- 1 )** | | 

DVR3 L | | MIXV,* | | ARG!2 I | (•)♦• I ? DVRBL | \ • < ITSPt - 1 XV • • || LHSA! 21 
I | M 1 /FI7SM1 • ,*1 *B, *0*B> I 

CALL CODE (OVRBL? ?TSS? | LHS A! 21 ? | •) = * | | PM | | 'SIGMA! I XV=2* ITS! * | I 
*! IXV-1 1 *• J | DVRBL? ?*IXV, * | | APGI21 | ?■>*•? | DVRBL I | * ! ITSPI-IXV, « 

| |LHSA(I 1 | j • 11/FITSMl • , • 1 • B . * O* B 1 1 
CALL CODE I DVR BL ? | TSS ? | LHS A (31 M *)-(* I | DVRBL If TSS I I LHSA! I) ? ? 

•) -SIGMA! I XV-2.ITS?* j | DVRBL? f » (I XV, • | ?LHSA!2 l| f • )** | | OVRBL? | 

*( ITSPI-IXV, • | ?LH$ A! 3 1 | |*))1/* | ? DVRBL? |*U,* | ?LHSA(2)| ? • )*, 

*1 *B • *0 * B 1 ; 

NMUL— NMUL+3 • NADD = NADD«-1 5 NMTS=NMTS+65 NMTL=MMTL-6; 

LZ! RETURN! *1 *B1 5 
END FUDGE? 


FUNCTIONS 
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* process; 

input; proc terror. cc>; 

/***************** ************* ***************** ************** ******* 

* PROCEDURE READS THE DEFINING SYSTEM CF DIFFERENTIAL EQUATIONS * 

* FROM THE SYSIN DATA SET * AND CHECKS TO SEE THAT THE EQUATIONS * 

* ARE BALANCED WITH RESPECT TO PARENTHESES A 

********************************************************************/ 

OCL 

BLNCD ENTRY T CHAR < * ) VARI RET UR NS T EIN FIXED)* 

COUNT ENTRY! CHAR T*> VAR.CHARTD) RETURNST8IW FIXED)* 

DELETE E NTRY < CHART *) VAR* CHART 1 ) ) * 

SPAN ENTRYT CHART *> VAR-* CHAP T * ) VAR* CHART * ) V AR* CHARI * ) VAR); 

DCL 

T E RROR »CC) BIN FIXED. LINE CHAR T 60 ) VAR. 

CWC2> CHAR <2 5) VAR I Nl T T * 01 FFERENTI ALEQU AT IONS • « 

•I NITIAL VALUES’) * CS CHART ADO ) VAR EXT* VS CHART 400) VAR EXT. 

I ED EXT. C48 CHART! ) EXT . T IVRBL.DVRBL ) CHART 4 ) VAR EXT! 

/* */ 
ON ENDFILET SYSIN) 

BE3IN; PUT EDIT!*** EOF READING SYSIN **• > TSKIPT2)*A>; 

GO TO lpa; 

END ; 

IF CC=1 

THEN DC; MRKR=0; 

PUT EDITT*** TAYLOR SERIES PROGRAM JAN. 1973* || 

• VERSION - LISTING OF INPUT EQUATIONS **•) TA); 

ERROR— 1 ; 

end; 

ELSE MRKR=I NDEXTCS. • A* t; 

put skip; 

lp: get edittline) tatso)); 

PUT EDITTLINE) T SKI P. COLUMN T4 ). AT 80) ) ; 

LINE- SUBSTR T L I NE. 1 ,72) ; 

CALL DELETETLINE, * •>; 

IF cc-.=© 

THEN DO; IF l NDEX TLINE.CWTCC) >-.-0 THEN GO TO LQ; 

PUT EDIT T * *♦ THE FOLLOWING CONTROL CARD IS INVALID ***• 
LINE) TSKIPT2) .A, SKIP, A); 


lpa : 

error=3; 

return; 


p 

o 

*4 

IF CC=I 



THEN DO; 

IF I NDEX TLI NE. ■ DP* >-»=0 
then do; ied=i; C48=*a*; end; 
else oo; ied=o; C48=»4*; end; 

IF INDEX TLINE.* T * )-=0 

then do; call sp ant l i ne. • t * ivrbl ) ; 

CALL SPANTLINE.* .•*•)• .DVRBLi; 

end; 


END 

4 * 


00000106 
00000200 
00000300 
00000400 
00000500 
00000600 
00000700 
00000800 
00000900 
00001000 
00001100 
00001200 
00001300 
00001400 
©0001500 
00001600 
00001700 
00061800 
00001900 
00002000 
00002100 
00002200 
00002300 
00002400 
00002500 
00002600 
00002700 
00002800 
00C 02900 
00003000 
00003100 
00003200 
00003300 
00003400 
00003500 
00003600 
00003700 
COO 03800 
00003900 
00004000 
00004100 
00004200 
00004300 
00004400 
00004500 
00004600 
00004700 
00004800 
00004900 


A- 19 



no : 


cc^o; put skip; go to lp; 

end; 

cs=cs f ) li ne; 

IF INDEX< LINE, • t* >=*0 THEN GO TO LP; 

put skip; 

cs-substr tcs • 1 , length ( cs)-nl 1 • ;* ; 

00 WHILE C WRKR<LENGTH«CS> ) ? 

MC = I NDEXISUBSTRI CS • MRKR + 1 > * * ; M * 
WS=SU0 STR < CS , MRKR+1 ,MC-1>5 


©ocosooo 

60005100 

©00032(50 

60005360 

©OCO540© 

OOOC5SOO 

©0005600 

COOO570© 

©0005800 


rort tv rt — it rv r*'r wv 9 

IF 8LNCOCWS5-*=© 

THEN DO! PUT EDIT<»** THE FOLLOWING EXPRESSION HAS AN 
•INCORRECT PAIRING OF PARENTHESES ***,WS) 
<SKIP<2) , A, SKIP, AI ; 

errors; 

end; 

IF COUNT C WS •■ = ■1-1=1 

THEN DO; PUT EDIT<*** SYNTAX ERROR IN THE FOLLOWING 
f | • **•, WS> ISKIPI25 • A, SKIP, A) ; 

E PROP =2 ; 

end; 

end ; 

END input; 


vvuvsvgo 
00606000 
I | ©00061©© 

00006200 
00606300 
0000646© 
00006560 
0000660© 
EXPRESSION *60006700 
©00C 6800 
©0006006 
0006700© 
©0007100 
60007200 
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* process; 

op: proc(Ch) returns<bin fixeoj; 

/ ********** 4 ****** *************************************************** 

* PROCEDURE CCWPARES THE INPUT CHARACTER CH WITH THE CHARACTERS * 

^ ♦» “* *» /» s » I* fi t» X * 

************************** *********************** *******************/ 
DCL CH CHAR < 1 > »I % OC <1 0 f CHARCII EXT? 

DO 1=1 TO io; 

IF CH=OC<I> THEN RETURN! 1 1 ; 

ENJi 

RETURN!©! ; 

END CP? 


©ooooioo 
©000020© 
©0000300 
©oe«©4©e 
00000500 
0000060© 
ooro© 7»© 

00300800 

0030090 © 

60001000 

00©011©6 

00001200 
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* process; 

optmze; proc!m,rmax> ; 

/* procedure eliminates redundant operations from the r matrix */ 

DCL 

8FDC ENTRY! BIN FIXED) RETURNS ! CHAR! I 5 ) VAR), 

IT *J,K,R,RMAX) BIN FIXED, IMP BIT! 1 1 ,M< *, * > CHAR! 15) VAR, 

I X ! 2 ,2 ) INIT!2,3 ,3,2) ; 

la: r-o; imp=*o»b; 

lab; p=r+i ; i =r+i ; 
lb; do whileo<=rmaxi ; 

IF M<R,l )-, = M!I ,1 ) THEN GO TO LD! 

IF M!R,1 )=*♦' I M!R, !)=•-* | M !R,1 )=•*• lM!R, !)=•/• 

THEN JMAX=2 ; ELSE JMAX=15 


ooccoirc 
00000200 
00000300 
90000490 
00000500 
00000600 
00 300 700 

ooocosoe 

©9000990 

ooooiqofl 

©ooonoe 

©0001200 

00001300 


LC 

LD 

LE 


/* 

rad; 


DO J=1 TC jmax; 

DO K=1 TO 2; 

IF M!R,K+1 >-=M!I,lX!K, J) > THEN GO TO LC, 

end; 

i mp— * i ' b; 

CALL RADII , P) " 

GO TO le; 
end ; 
i = i+i ; 
end; 

IF R<RMAX— T THEN GO TO LAB; 

IF tMP THEN GO TO LA; 

PROC! I ,J> ; 


©990249© 

0000150© 

ooeci 600 

©9001700 

00001600 

©0001900 

00002000 

0000210© 

©006220© 

00002300 

©0002400 

00002500 

4/00002600 

00002700 


/ ** ******* ************************************** ******* *** * *********** 0000280© 


* PROCEDURE DELETES RCW I, AND REPLACES REFERENCES TO ROW I WITH * 0000290© 

* ROW JIN THE RECURRENCE MATRIX * ©0003000 


********************************************************************* /©0003100 


DCL I ,J,K,L,RCW BIN FIXED! 

K=I+1 ; 

DO WHILE! K<=R MAX) ; 

M( K-l ,1 )=M!K,l ) ; 

DO L=2 TO 4; 

IF SUBSTRI M ! K *L ) ,1 ,1 * ? • 

then do; m!k-i,li=mik,l); go to lf; end; 

RGW=SUBSTR(M!K,L) *21 ; 

IF ROW=J THEN ROW=j; ELSE IF ROW>! THEM ROW=ROW-1; 
M ( K— I , L ) = • 7 • (I BFDCCROW); 

lf; end; 


k=k+i ; 

END ; 

RMA X=RMAX-1 ; 
END rad; 
end optmze; 


00003200 

00063300 

00003400 

60003500 

00003690 

00603700 

©0003890 

00003900 

00004090 

00004100 

C0004£©0 

00004300 

C0004400 

00004590 

00004690 

00004790* 
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* process; oooooioo 

rmat: proc<cs.r.rmax.ko> recursive; 00000200 

/ ******************************************************************** 00000300 

* FACTORIZATION OF THE DIFFERENTIAL SYSTEM INTO CANONICAL FORM * 00000406 

* USING THE ALGORITHM DESCRIBED IN « BOUNDED CONTEXT TRANSLATION * . * ©0000500 

* BY R. GRAHAM* AFIPS-SJCC V 25 (19641* P. 21 4 00000600 

******************************************************************44/ 00000700 


/* 


LI : 


L2: 

L3: 


DCL 

OP 5 NTRY ( CHAR ( 1 ) ) RETURNS (BIN FIXED). 

BFDC ENTRY! BIN FIXED) RETURNS C CHAR( 15 ) VAR )• 

REPLACE ENTR YC CHAR ( * ) VAR.CHAR<4> VAR.CHAR(4» VAR.BITCll). 

SFNL ENTRY(CHAR(4) VARI RETURNS(BIN FIXED). 

SPAN ENTRY(CHAR< 4) VAR.CHAR(4) VAR.CHAR(4) VAR.CHARI4) VAR I ? 

DCL PCS!) BIN FIXED STATIC IN IT< 6. 6. 7, 7, 5. 3. 4 • 1 • 2. 8. 8) . 

SC 500) CHAR < 1 51 VAR INTT ( < 500) L( 250) CHAR ( 1 5 ) VAR* 

TYPE BIN FI XED.R ( 4*4) CHAR<1S) VAR.CS CHAR <4 ) VAR. 

WT C H AR( 400) VAR .FN CHAPC4) V AR.FTBLC 100. 2 ) BIN FIXED EXT. 

WS CHAR(400) VAR .R MAX BIN FIXEO.DVRBL CHAR <4 ) VAR EXT! 

DCL CST(ieO) CHAR (25) VAR EXT.IC EXT.KO 9IT(1).IED EXT. 

NEST BIN FIXFO STATIC INlf (f» >, DEBUG BITCH EXT.NEO EXT. 

PMD CHAR (36) VAR FXT? 

4/ 

ko= u * b; 

NEST=^ST-fl ; 

IF NEST=1 C NEG=1 C DEBUG 

then do; 

put page; 

PUT EDIT (* RMAT ENTRY LEVEL K TYPE (C. O. V | E)«* 

* R OP*.«A(l )«.* A(2 )• .*A(3I» I <SKIP(2).A.SKIP<1 ),A 
XC33) .A.XC5) *A.3 (X(ll).A)); 

END ; 

DO 1=1 TO LENGTH! CS 1/45+1 ; 

PUT EDIT! SUBSTRCCS.l ♦ C I -1 ) 445. M IN C 45.LENGI H( CS )-< I - 1 ) 445 ) ) ) 

( SKIP.XCl 5) .A) : 

END ; 

j=i ; 

k=2 ; 

MRKR=2 ; 

S(1 )=SUBSTR(CS.l .1) ; 

type=8; 
lci ) = su ) ; 

IF SC K) = • • THEN CALL CHECK? 

IF DEBUG THEN PUT EDI T C NEST. K. TYPE *S C K ) ) (SKIP. 3 FC4). X( 3 ) * A ) « 

IF TYPE<I THEN GO TO L25 
IF TY PE“*=6 THEN GO TO LAI 
j=j+1 ; 

LC J) = S(K) ; 

k=k+i ; 

IF K>DI MC S . 1 ) 


00000800 

0OOCC9O0 

ooootooo 
00001100 
90001200 
00001300 
00001400 
OOCOISOO 
00001600 
00001700 
00001 800 
00001 900 
00002000 
00002100 
OOOO220C 
00062300 
00002406 
00002506 
000 62600 
000C2706 
00002800 
00002906 
00063000 
00003100 
00003206 
C060330C 
00003400 
00603500 
00003606 
00063700 
00003800 
00063906 
00064006 
00004100 
00064200 
600C4306 
00004400 
00004506 
00064600 
00004700 
60064806 
00064906 
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THEN DO; PUT EDIT<»** OVERFLOW IN S TABLE **• * <SKIP(2),A); 

ko»’Q*b; return; 

end; 
go to Li ; 

la: IF P<OP<L<J-lin<P<TYPEl THEN GO TO L8; 

RMA X=RMAX*1 ; 

IF LC J—l 1 * THEN GO TO LAB; 

IF LENGTH<L< J-2) > <LENGTH<DVRBLI 1 SUBSTR <L< J- 2 J* 1 • LENGTH! OVRBL } ) 
-t=DVRBL THEN GO TO LAA; 

CALL SPAN CL? J-?) . »• .FN) * 

R(RMAX«1 > = •*» ; 

r<rmax*2) *dvrblM • < • Hen 1 1 *»• ; 

R(RMAX*3> -L< Ji ; 

R<RMAX.4I=R<RNAX*2> ; 

GO TO lac; 

laa: IF SUBSTR <L(J)*1*1 )-*=’ ?• | SUBSTRCLC J »*2J-*=BFOC<RMAX-l > 

THEN GO TO LAB; 

RMAX=RMAX-1 ; 

RCRMA X,4) -L< J-2) ; 

GO TO LAC ; 

lab: IF U = G I NDEX <L < J ) * * A* )-*”0 

THEN DO? IF CST<SUBSTR<L< J>*2 ) 2* THEN GO TO LAZ; 

R< RM AX » 1 » = »*'; 

R<RMAX*2)— L( J— 2 } * 

R<RMAX*3)=L< J-2 I J 

end; 

ELSE 

laz; do; r<rmax*i »=l< j-i i; 

P<RMAX *2 )=L< J— 2> » 

filDUiV .3 > -l_ / J * • 


©oocsoe© 

©0005100 
0000 5200 
00005300 
00005400 
00005500 
COO 05600 
CO0057OC 
00005800 
nnnr 5 on it 

00006000 

09006100 

00006200 

00006300 

00006400 

00006500 

00006600 

00006700 

00006800 

00006900 

00007000 

C00Q7100 

00007200 

00007300 

00007400 

©0007500 

00007600 

©0007700 

00007800 


end; 

IF L< J-l ) = •=* 

THEN R<RMAX*4 )=R< RMAX*2) ; ELSE R< RMAX * 4 1= • 7 • | 1 BFDC< RM A X ) ; 
LAC: IF DEBUG THEN PUT EDI T < RMAX* < R< RMAX* M ) DO M=1 TO 4)> 
<SKIP<X<60I ,F<2> ,X(1 ) «A(2)*X <5 1*4 A<15»»; 
j= j-2 ; 

L< J)~ * 7 • || BFOC<RMAX>; 

GO TO la: 

lb: if TYPE=7 THEN DO * L(J-ll=L<Ji; J= J-l ; GO TO L3* END* 

IF TY PE-»-9 THEN 60 TO L2 : 

IF K- 2 THEN CS=S(2»; ELSE CS= • ? • | | BFDC< RMAX I ; 

NEST=NEST-1 ; 

/* 

check: proc recursive; 


00008000 
00008100 
00008200 
00008300 
00008400 
©000050© 
0000860© 
00008700 
00008800 
00008900 
09099000 
00009100 
*/ 00009200 
©0009300 


/******************************************************************** 00009400 

* PROCEDURE SCANS A SEQUENCE OF CHARACTERS BEGINNING WITH MRKR TO ♦ ©0009500 

* DETERMINE WHETHER THEY SPECIFY A CONSTANT* VARIABLE* OR OPERATOR * 000096©© 
********************************************************************/ ©900970© 


DCL 


©0099800 
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3LNCD ENTRY < CHAR C * ) VAR) RETURNS (BIN FIXED)* 

FLIP BIT (II* CH CHAR ( 1 ) *STRI NG CHARC15I VAR INITC**)* 

BREAKS ENTRY<CHARC*I VAR*CHAR<*> VAR*CHAR< * ) VAR)* 

BREA KF E NTRV < CHAP C *) VAR*CHARC*I VAR,CHAPC*> VAR); 

/* 

CH=SUBSTR (CS* VRKR,1 )? 

sc K)=Ch; 

MRKRs MRKR-H ; 

IO=DPCCH> ; 

IF I 0—0 THEN GO TO LG? 

IF < T YPE — 5 | T Y PE=6 1 & I O < 3 

THEN IF VERIFVCSUBSTR<CS*MRKR* 1 )«• *9123456769* )-*=0 
THEN GO TO LGA » ELSF GO TO LJI 
IF TYPE >0 CTYPE-=6 &T YPE-»=7& IO-»=6&IO-»=7 
THEN GO TC LJ; 

IF 10=9 THEN DO? TYPE=IO? GO TO NO? END • 

ch=substrccs*wpkr*i i ; 

NOP=OP CCH ) ; 

IF NOP=3 & 10=3 THEN DOS S<K)=S<K) \ | CH? MRKR=MRKR4’l 5 END* 

IF LE NOTH C S ( K 1 ) —Z THEN TYPE =115 ELSE TYPE =105 
GO TO nd; 

LG: 1= VERIFYCSCK) **ABCOEFGHIJKLMNOPORSTUVWXYZ?« l-»=0 THEN GO TO LJ5 

/* PROCESS A VARIABLE */ 

lga: fhp=»9'8; type^o; »5=$<k>; 

Lh; ch=substr ccs * mrkr *1 1 ; 


©0009909 
00010006 
00010100 
09010200 
*/ 00010300 
00010400 
C0010500 
09010000 
00010700 
00010890 
00010900 
00011000 
00011100 
00011200 
00011300 
00011400 
00011590 
00911600 
00011700 
09011800 
00911900 
00012000 
09012100 
90012290 
00012390 


NOP=OP CCH I * 

IF N0P=0 THEN GO TO LI 5 

IF NOP-t-6 t N C P-*= 7 6 -.FLIP THEN GO TO LIA5 
IF BLNCDC VS | fCHXO THEN GO TO LIA5 
FLI P= • I »B5 
LI 3 *S=WS | f CH 5 

MRKR=MRKR+1 * 

IF -FLIP J BLNCDC WS )-*=9 THEN GO TO LH* 

CALL SPANCM* l|WS.')»*M , ,WT); 

IF WT=0 VRBL THEN 

lia; do; sck)=»s; go to nd; end; 

lfn=sfnl< vs) ; 

IF LFN=0 THEN GO TO LIA; 

CALL 9REAKFf*S*M • * WT | ; 

CALL SREAKBCWS* •) • *NT ) ; 

IF VERIFY CVS* PMD) =0 THEN GO TO LIB* 
WS=**TEMPN=* \ iMSll * ** 5 
CALL RMAT fi*RMAX *KQ) ; 

IF-KO THEN RETURN; 

PCR MAX*4) =•?• i IBFOCIRMAX) : 
lib: CALL SFNCILFN*KS.R*RNAX*S< K),FTBL*K015 
IF -»KO then return; 

GO TO no; 

PROCESS A CONSTANT */ 


09012490 
00012590 
00C 12600 
00012790 
09012800 
00012990 
00013090 
00013100 
©0013200 
00013300 
00013400 
00013500 
00013600 
©0013700 
00013800 
00013900 
00014090 
©0014100 
0001420Q 
00014300 
00014400 
00014500 
00014600 
00014700 
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LJt TYPE *— 1 * 0001480© 

IF IC^IO© 00014900 

THEN DOS PUT EDIT !*♦* OVERFLOW IN CONSTANT TABLE CSK IPt 2 J *A) ;©0CIS0©C 

KG=*©*e; return; oooisio© 


end; 
ic=ic+i ; 
cstuc>=s<kj ; 

lja: ch- substr (Cs«mrkr-»ii; 

IF VERIFYICH, ‘.©l 23456789E04— * >-.=© THEN GO TO LK; 

I F ? CH= * }CM=«-« ) &VER IFY ( SUBSTR < CST < IC I* LENGTH! C5TCIC))«1)« 
•ED • ) -,-© THEN GQ TO LKJ 

IF (CH*»E* | CH^*D*> S VERIFY {STRING, * *1234567890* >-*=© 

THEN GO TO LK*, 

C ST ( I C > — C ST (IC) If CH I 
STRING-STRI NG ( | CH ; 

mrkr=mpkr*i ; 

GO TO lja; 
lk: if i c -l 

THEN DO; if lED^l & INDEX(CST (IC).'.' l-=© 

THEN do; call replaceicst UC>* 'E»* *0* 

IF I NOEX <CST (IOt'D' 

then cstiici=cst<icM I* do»; 

end; 

do i =1 to ic-i ; 

IF CST <11 -»=CST I IC I THEN GO TO L J3 * 

s<k> = * \ | 8foc<ii ; 

I c-tc-i : 

return; 

ljb: end; 

END * 

S(K )=•*• tf BFOCIIC); 

return; 

nd: end check; 

END rmat; 


©0015200 
©0015300 
©0015400 
©0015500 
© 001560 © 
00015700 
©0015800 
©0015900 
©0016000 
©0016100 
00016200 
00016300 
©0016400 
00016500 
©0016600 
00016700 
©0016800 
00016900 
©0017000 
00017100 
©0017200 
©0017300 
00017400 
00017500 
00017600 
©00 17700 
©0017800 
©0017900 

cooieoco 

©0018100 
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* process; 

SFNC: PROC 1LFN.WS* R.RMAX *S * FTBL.KQ > ; 

/******************************************************************** 

* PROCEDURE INSERTS THE PROPER ENTRIES IN THE R MATRIX FOR AN * 

♦ ALLOWABLE FUNCTION REFERENCE * 


00004100 

00000204 

04000300 

00000400 

eoooosoo 


** ** ******* ****** ****** *** ******* ************** ********** ***********/ 00000600 


/* 


LA 


LB 


END 


DCL 9 F 0 C ENTRY! BI N FIXE 01 RETURNS! CHAR ! 151 VAR I.RMAX BIN FIXED* 
R(*.*> CHAR f *> VAR*S CHAR! * ) VAR.WS CHAR! * I VAR* 

CLFN CHAR ! 15 I VAR.FTBL C*. *> BIN F IXED.KFMAX EXT* DEBUG BIT! 1 > E 
KFL ! 10 1 INIT! 1 , 2 . 3 * 4 * 4 . 4 . 7 . 7 . 7 . 1 © ). 

K-U! 10 ) I NIT ! 1 . 2 , 3 * 6 ,6 * 6 , 9 . 9 , 9 ,10 II 

CLFN=BFDC!LFN) ; 

DO 1=1 TO RMAX; 

IF R! 1 * 1 ) -*= * X • j R ! 1 . 2 ) -*=CLFN f R!t* 3 )~.=WS 
THEN GO TC LAS 

ELSE DO? S= •? * f | BFOC(I); RETURN; eno; 

END ; 

DO K= KFL! LFN| TO KFUfLFNl; 

rmax=rvax+ i ; 

IF K =10 

THEN DC; RIRMAX.l )= ** 4 * ; 

R!RMAX* 2 I=XS; 

R ! RMAX * 3 ) = • #1 * ; 

GC TO LB; 

end ; 

Rf R MA X * 1 )=•*•; 

R< RWAX* 2 )=BFDC!K) ; 

R( RMAX * 3 ) = MS; 

R! RMAX * 4 ) = • ? * If BFOCfRMAX); 

IF DEBUG 

THEN PUT EDIT! RMAX « < R! RMAX * Ml DC M =1 TO 4)1 ( SK IP. X < 60 ) * F < 2 ) 
X!l) ,A (21 *X!S) *4 A( 15 l); 

IF K=LFN THEN S=*?* | f eFDC! RMAX) ; 

END ; 

KFMAX=KFMAX +1 ; 

IF KFMAX>DI M(FTBL* 1 ) 

THEN do; PUT EOIT!'** OVERFLOW IN FUNCTION TABLE **• ) 

I SKI P! 2 ) * A) ; 

kc=»o*b; return; 

end; 

FTBL! KFMAX* 2 ) =RMAX; 

FTBL! KFMAX.l ) =RMAX-KFU !LFN ) fKFL <LFN) ; 

SFNC ; 


©0000700 

eooooaoo 

XT *00000900 
00001000 
ooeonoc 

*/ 00001200 
00001300 
00001400 
00001500 
00001600 
00001700 
000 01 800 
00001900 
00002000 
00002100 
00002200 
00002300 
00002400 
00002500 
00002600 
00002700 
00002800 
00002900 
00003000 
00003100 

. 00003200 

00003300 
00003400 
00003500 
00003600 
00003700 
04003800 
000C390C 
00004000 
00004100 
00004200 
00004300 
00004400 
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* process; 

SFNL: PROCIWS) RETURNS (BIN FIXED!? 

/* PROCEDURE OETER WINES WHETHER W$ IS AN ALLOWABLE FUNCTION */ 

DCL 

WS CHARI *1 V AR»NF EXT *F CHAR <6 I VAR*FNC16*2» CHAR <6! VAR EXT; 
FssSUB STRC WS*1 *1 NDEX <WS * ■ t • 1-11* 

oo i* i to nf: do u-i to a; 

IF F=FMI,J! THEN RETURN! t 1 * 

end; eno$ 

RETURN' A! : 

ENO SFNL? 


©0*00100 

0060026 * 

00006300 

*6400400 

©0*065©* 

60000600 

000007*0 

006008*6 

00000960 

60001000 

66C01100 
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* process: 

SIGMA: PROCfSTRINGl RECURSIVE: 

/* PROCEDURE CONSTRUCTS A FORTRAN DO LOOP FOR A SUMMATION OPERATOR 
DCL 

B=DC ENTRY! 0 IN FIXED I RETURNS ! CHAR! 1 5 ) VAR I* 

8REAKF ENTRY (CHAR!*! VAR* CHAR (11* CHAR! * ) VAR »* 

CODE ENTRY<CHAR<*> VAR « 81 T ( l 1 , SIT < 1 ) ) . 

EXTRACT ENTRY(CHAR!*I VAR*CHAR(*I VAR.CHAR!*! VAR»* 

DCL EXA CHAR!*©©) VAR • NSGM A 0IN FIXED STATIC EXT, 


000<J01f>© 
00000200 
*/ 00000300 
00000400 
00000500 
00000600 
00000700 
00000800 
00000900 


(STN*SMN) CH AR(5 ! VAR,EIS(3I CHAR(8> VAR EXT « STR ING CHAR!*! VARI OOCOIOOO 


/* 

NSGM*=NSG*A+1 : 


*/ 00001100 
00001200 


SMN"=8FDC! NSGM Al : 


00001300 


STN=8FDCMOOO 4NSGMAI ? 


00001400 


CALL BREAKF (STRING* • : * .EXAM 00001500 

CALL CCDEC*SGMA*1 | SMN || *=0 .0* * • 0 • B ■ « © » B 1 1 00001600 

CALL CODE ( *0O * | j STM | | • • f | EX A* • 0 * 0* • 0 • 81 : 00001700 

CALL CODE ( STN f ( •SGMA* f (SMNf | » =SGMA* j | SMNf f « +! « f | SUSSTR ( STR ING * I * 00001800 
LENGTH! STRING! -I ) H * > * * « 1 • B* ■ 0 • 0> I 00001900 


END SIGMA; 


00002000 
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* process; ©eaooioo 

stint: proc; oo©«©2©« 

/* PROCEDURE INITIALIZES ALL EXTERNAL BIT AND CHARACTER STRING VBL • S */©0000300 
DCL < 00000400 

NF I NITC 10) • CIVR8L*DVflBL) CHAR (4) VAR.EISC3) CHAR <81 VAR* 00000500 

OCIIOI CHAR* 11 *FNC13*2) CHAR C6 > VAR*SVSA BITCIKPMD CHARC36) VAR 00000600 

) ext; 0000070C 

/* */0(M5Cf80<! 


sysa=*i *b; 

PM!5 = • ABCDEFGH I JKL NNOPORSTUVWX V ZI 234567890 * S 
I VR DL — * T ~ * uV wot— • ¥•' i 


EISU )=*EQUATICN* ; 

EIS <2 ) = • INTEGRAL* ; E IS C 3 >= * SI GMA* ; 

OCC 

i >=•♦ • ; 

OCC 2) 

— i.i 

; occ 3)=***; occ 4)=*/«; occ 5)=*=*; 

OCC 

6>=*c • ; 

OCC 7) 

= *) * 

; occ 8)=***; occ 9)=*s*; occto>=*%*; 

FN( 

i ,1 > = 

•Exp* ; 

FNC 

1.21= • OEXP* ; 

FNC 

2*1) = * ALOG1 0 * * 

FNC 

2*2 > = • DLOGIO* ; 

FNC 

3*1) — 

* ALGG* ; 

FNC 

3*21= * OLOG* ; 

FN< 

4*11 = 

• si n • ; 

FNC 

4*2)= * DSIN* 5 

FN< 

5 *1 ) = 

•cos* : 

FNC 

5,2)= *dcos*; 

FNC 

6*1) = 

•TAN* ; 

FNC 

6*2)= * DT AN* « 

FN( 

7*1 > = 

• S I NH • : 

FNC 

7,2)= *OSINH*; 

FN( 

8*1 > = 

•CCSM* ; 

FNC 

8,2)= *OCOSH»? 

FNC 

9*1 > = 

•TANH* ; 

FNC 

9,2)= *otanh*; 

FNC 10 *1 ) = 

•SORT* ; 

FNC 10 *2 ) = * DSQRT* 5 


END stint; 


0000090© 
000010©© 
00001100 
00001 200 
00001300 
00001400 
00001 500 
00001600 
00001700 
00001 000 
00001 ROC 
00002000 
©0002100 
00002200 
00002300 
©0002400 

©oooesoo 
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♦ process; 

BFOC5 PROC<N) RETURNS (CHAR! 15) VAR)* 

/************ *4444444*4*44**4444*4*444****4*444*444 ****************** 

4 BIN FIXED NUMBER N IS CONVERTED TO CHARACTER STRING WITH ALL * 

* BLANKS RESULTING FROM THE CONVERSION DELETED * 

************** ******* ************ ***************************44*44 4*4/ 

DCL DELETE ENTRY ( CHAR ( *1 V AR • CHAR ( III * CS CHAR<15) VART 

CS=CHAR(N) 5 

CALL DELETEtCS* • ')? 

RETURN(CS) ; 

end bfdc; 


00000100 

00000200 

00000300 

00000400 

00000500 

00000600 

00000700 

00000800 

00000900 

00001000 

00001100 


♦ process; 

BFTC: P ROC ( BF * I ED ) RETURNS ( CHAR ( 50 ) VAR); 
/******************************************************************** 

4 BIN FLOAT ( 53 ) NUMBER BF IS CONVERTED TO CHARACTER STRING WITH 4 
4 ALL BLANKS RESULTING FROM THE CONVERSION DELETED * 

********************************************************************/ 
DCL 

DELETE E N TRY C CHAR (4) VAR .CHAR ( 1 >> , 

REPLACE ENTR Y(CHAR<4) VAR.CHAR<4> VAR,CHAR(4) VAR,BIT<1)1, 

BF BIN FLOAT ( 53 ) »CH CHAR ( 50 ) VAR.SF BIN FLOAT * 

IF IE0=0 THEN DO; SF— BF ; CH=SF; END* ELSE CH=BF» 

CALL DELETE! CH» • » ) ; 

IF IED=1 THEN CALL REPLACE!CH • • E « , » D« • »0 ■ B ) ; 

RETURN(CH > ; 

END BFTC ; 


00000100 
00000200 
00000300 
00 000400 
00000500 
00000600 
00000700 
00000800 
00000900 
00001000 
00001100 
00001200 
00001300 
00001400 
00001500 


* process; 

blncd: procccs) returnscbin fixed); 

/****************** ***************************** 44 * * 4 ************* *** 

4 THE DIFFERENCE BETWEEN THE NUMBER OF RIGHT ANO LEFT PARENTHESES 4 
4 IN THE CHARACTER STRING CS I S RETURNED 4 

*************************************44444*44444444%****************/ 


DCL CS CHAR! 4 ) VAR, C CHAR(l) VAR, <IL*IF) INtT(O), I? 
DQ 1=1 TO LENGTH (CS) ; 

C=SUBSTR(CS,I ,1 ) ; 

IF C=»(» THEN IL=IL*15 
IF C=*)« THEN IR=IR+i; 


END ; 

RETURN (IL-IR); 

end blncd; 


00000100 
00000200 
00000300 
00000400 
00000500 
00000600 
00000700 
00000800 
00000900 
00001000 
000011 00 
00001200 
00001300 
00001400 
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* process; 

breaks; procccs »c*bs> ; 

/**********♦********************************************************* 

* THE CHARACTER STRING TO THE RIGHT OF THE CHARACTER VARIABLE C IN * 

* THE CHARACTER STRING CS IS PLACED IN THE CHARACTER STRING BS * 

* C || BS IS DELETED FROM CS * 

* example: cs = »abcde' * 

* CALL BREAKS CCS • 1 C * *BS ) * 

* GENERATES C S - * AB * • BS = * DE * * 

****** A*** ******************************************************* ***/ 

DCL CS CHARC*) VAR* C CHAR(1>. SS CHAR C* ) VAR* IX? 

bs= * * ; 

DO I “LENGTHCCS > TO 1 BY -I ; 

IF SUBSTRCCS. I • 1 > -C THEN GO TO LA? 

IF I+1CLENGTHCCS) THEN BS=SUBST R < CS • I * 1 > * 

IF I>1 THEN CS=SUBSTR(CS, l *1-1 ) ; ELSE CS=»»J 

GO TO lb; 
la; end; 

LB : END BREAKB : 


00000100 
00000200 
00000300 
00000400 
00000500 
00000600 
oooootoo 
00000800 
00000900 
00001000 
00001 1 00 
00001200 
00001300 
00001400 
00001500 
00001600 
00001700 

ooooieoo 

00001900 


* process; 

breakf: procccs *c* BS) ; 

/*** ********* ********** ** * * *** **************** ******* ************* *** 
* THE character STRING TO THE LEFT OF THE CHARACTER VARIABLE C IN * 


* THE CHARACTER STRING CS IS PLACEO IN THE CHARACTER STRING BS * 

* BS || C IS DELETED FROM CS * 

* example: cs = ■ abode* * 

* CALL BREAKFCCS* *C< *BS) * 

* GENERATES CS = »DE<* BS = * A8 • * 


***** ********* ********************************************* *********/ 
DCL CS CHARC*) VAR. C CHAR(l). BS CHARC*) VAR.IX; 

bs=* * ; 

I X = I NO EX(CS»C)-t ; 

IF IX=-1 THEN RETURN; 

IF IX<1 THEN GO TO LA? 

BS=SUBSTR<CS. 1 .1 X) ; 

LA: if I X+2<=LENGTH( CS ) THEN CS=SUBSTR C CS * IX + 2 ) ♦ ELSE CS=**; 

LB : END BREAKF ; 


00000100 
00000200 
00000300 
00000400 
00000500 
00000600 
00000700 
00000800 
00000900 
00001000 
00001 ! 00 
00001200 
00001300 
00001400 
00001500 
00001600 
00001700 
00001800 
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* process; 

count: procics.C) returnsibin fixed); 

/* ****** ********** *************************************************** 

* THE NUMBER OF TIMES THE CHARACTER STRING C APPEARS IN THE * 

♦ CHARACTER STRING CS IS RETURNED * 

*** ******* ^♦♦♦f^t************************************************ AAA/ 

DCL CS CHARI A ) VAR* C CHAR I * ) VAR . 1 MRKR. I C ) INITIO); 

la: ixc= inoexisusstr(cs*mrkr+i ) *c> ; 

MRKR=MRKR+ I XC ? 

IF IXC-*=0 THEN IC-IC + i; 

IF IXC=0 | MR KR= LENGTH ( CS ) THEN RETURN! ICM 

go to la; 
end count; 


00000100 

00000200 

00000300 

00000400 

00000500 

00000600 

00000700 

00000800 

00000900 

00001000 

00001100 

00001200 

00001300 


a process; 
delete: procics.o; 

/A A A A A AAAA A A**’* A A A A AAA A AAA A A A* A AA A A AAA A AAA AAA A AAA A AAA A AAA A AAA AAAAA AAA 
A EXAMPLE: CS = 'AOCDE* CS = * A * * BCD * * E * A 

A CALL DELETEICS * *C • ) CALL DEL ETE I CS • * * • • I A 

A GENERATES CS * • A8DE • CS = « AE • A 

********* ******** * ************ AAA A AAAA AAAA AA AAAA AAA A AAA* AAAA AAAA A ***/ 
.DCL C CHAR II)* CS CHAR I A ) VAP . I INITI2). 

FLIP BIN FIXED I N I TC 1 ) .SY M CHARI 1 ) INITI ,,,, i; 

CS=* •|| CS || * * ; 

IF C-SYM THEN GO TO LH ; 

DO WHILE! I <LENGTH<CS) ) ; 

IF SUBSTR (CS.I.l) =C 

THEN CS=SUBSTRICS. 1 .1-1) || SUBST R ICS ♦ 1+ 1 ) ; 

ELSE IslAi; 

end; 

go to li a; 

lh: DO WHILE! I <LENGTH(CS) ) ; 

IF SU8STRCCS.I .1 >=SYM 

then do; flip=-flip; GO to lha; end; 

IF FLIP<0 THEN 

lha: do; cs-substrics, i * 1 - 1 > II substrics.iai >; go to li; end; 

ELSE 1=1*1 ; 

L I : end ; 

LI A : CS = SUBSTR( CS . 2 .LENGTH ( CS ) -2 ) l 
end delete; 


00000100 

00000200 

00000300 

00000400 

00000500 

00000600 

00000700 

oooooeoo 

00000900 

00001000 

00001100 

00001200 

00001300 

00001400 

00001500 

00001600 

00001700 

00001800 

00001900 

00002000 

00002100 

00002200 

00002300 

00002400 

00002500 

00002600 
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* process; 
libf: proc<cs* ied) ; 

/*** 4c* 4c 4 ^ 4t***# ***************************** ****************** ******** 

* IF I ED = 0 ALL DOUBLE PRECISION FORTRAN LIBRARY FUNCTIONS IN THE * 

* CHARACTER STRING CS ARE REPLACED WITH SINGLE PRECISION FUNCTIONS * 

* AND VICE VERSA IF IED =1 * 

4c ***** ****** * * **4c* ********************* ★ * * * ** ***** ********* *** 4+ * / 

DCL IED, CS CHAR<*> VAR ,NF INITI25), 

REPLACE ENTRY (CHAR < *) VAR,CHAR<*> VAR.CHAR<*> VAR*SIT(D); 

DCL FN<25,2> CHAR ( 6 ) VAR INIT( 


'EXP* • 

•DEXP* , 

* ALOG 10 • • 

• DLOGIO •* 

* ARSIN* , 

• DARS IN • « 

• ARCOS • * 

•DARCOS* , 

•AT AN* • 

•DAT AN* , 

* AT AN2 • » 

• DAT AN 2 * , 

•SIN* , 

* DSI N • • 

•COS* * 

• DCOS • , 

•TAN* , 

* DT AN • • 

•COTAN • * 

•DCOTAN • • 

•SORT • • 

•DSORT •• 

• TANH* • 

•DTANH* , 

• S I NH • , 

• DSI NH • * 

•COSH* • 

• DCOSH • * 

•ERF • , 

• DEPF* t 

•ERFC * • 

♦ DERFC • • 

•GAMMA* , 

•OGAMMA* , 

* ALGAMA* , 

* DLGAMA • , 

•AMOD • • 

•DMOO • , 

• ABS* • 

•DABS •* 

• AM AX 1 • • 

•OMAX1 • • 

•AMINl • , 

•OMIN1 *, 

•FLOAT • , 

• DFLOAT • • 

• S I GN • , 

* DS I GN • , 

• ALOG • • 

* OLOG • 

> ; 





DO 1=1 TO nf; 

CALL PE PLACE! CS*FN( I • 2-1 ED ) ,FN( I , 1 4- IED > , • 1 « B) ; 
END ; 

END LIBF; 


oooqoioo 
00000200 
00000300 
00000400 
00000500 
00000600 
00000700 
00000800 
00000900 
00001000 
00001100 
00001200 
00001300 
0000 1400 
00001500 
00001600 
00001700 

poooieoo 

00001900 

00002000 

00002100 

00002200 

00002300 
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* process; 

replace: prqc ! c s* sa * sb * opi ; 

/***♦************* ******** ******************************* ******** * *** 
* ALL APPEARANCES OF THE STRING SA IN THE CHARACTER STRING CS ARE * 


* REPLACED WITH THE STRING SB IF THE BIT(l) VARIABLE OP = *0*0 * 

* IF OP - * 1 * B ONLY THOSE OCCURENCES OF SA SOUNDED ON THE RIGHT BY * 

* THE NULL CHARACTER OR ♦-/*>! ARE REPLACEO WITH SB * 

* example: CS = *A*BA» * 

* CALL REPLACE (CS ♦ *A * * * Z * * *0*81 * 

* GENERATES CS = *Z«-BZ» * 

* WHILE CALL REPLACE (CS * * A * , • Z * * * l ' B ) * 

* GENERATES CS = f Z+BA« * 

* restriction: the CHARACTER X MAY not appear in CS* SA* OR SB * 


*************+******♦*+********* * * *************************** ***** *** / 

DCL 

CHECK ENTRY (BIN FIXED) RETURNS! B IT < 1 )) * 

VERIFY ENTRY<CHAR{*> VAP.CHAM*) VAR) RETURNS (BIN FIXED)* 

{ CS • SA * SB ) CHAR ( * ) VAR*OP BITtD.OR BIT(l) INIT<*0»B); 

/* */ 
LS A- LENGTH ( SA ) 5 

la: mrksa=index<cs*sa> ; 

IF MRKSA=0 THEN GO TO L8 ; 

IF -»QR THEN DO * CS=* •JJCSJJ* M MRKSA=MRKSA4- 1 ? END* 

Q R — * 1 * B : 

CS=SUBSTR(CS.l .MRKSA-1 ) || * %' || SUBSTR ( CS * MRKSA+LS A ) ; 

GO TO LA I 

LB: IF -iQR THEN RETURN? 

mrksa=inoex<cs* •%*) ; 

IF MRKSA=0 THEN GO TO LC ? 

IF -.OP | CHECK(MRKSA) 

THEN CS=SUBSTR<CS.l .MRKSA-1 ) || SB || SU8ST R ( CS * MRKSA+1 ) ? 

ELSE CS=SU9STR(CS. 1 *MRKSA-1 > j| SA jj SUBSTR C CS . MRK SA «■! ) J 
GO TO lb; 

LC : C S — SUB STR ( CS » 2 .LENGTH! CS ) —2) ? 

/* */ 


00000100 

00000200 

00000300 

00000400 

00000500 

00000600 

00000700 

00000800 

00000900 

00001000 

00001100 

00001200 

00001300 

00001400 

00001500 

00001600 

00001700 

00001800 

00001900 

00002000 

00002100 

00002200 

00002300 

00002400 

00002500 

00002600 

00002700 

00002800 

00002900 

00003000 

00003100 

00003200 

00003300 

00003400 

00003500 


check: proc ( ix) returns!BIT< i ) ) ; 

DCL IX BIN FI XED. ( IL.I R) BIN FIXED I NT INIT<0>; 

IF IX-1>0 THEN IL = VERIFY<SU8STR(CS*IX-1. 1 ) , «=(♦-*/ )« ) ; 
IF IX*1<=LENGTH(CS) 

THEN I R= VERIFY! SUBST R ( CS » I X+ 1 , 1 ) . *=<+-*/ ) • ) ? 

RETURN ( IL+IR=:0) ? 

END check; 

END replace; 


00003600 

00003700 

00003800 

00003900 

00004000 

00004100 

00004200 

00004300 
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* process; 

span: p«oc<cs.sa#sb.sc> ; 

/**%********* ************************************************** ****** 

* THE CHARACTER STRING IN CS SPANNED BY SA AND S0 IS RETURNED IN SC* 

* EXAMPLES CS = «ABCDEF» * 

* CALL SPAN! CS* * AB * * *F • »SC ) * 

* GENERATES SC “ ■CDE* * 

************** ******************* ******************************** ***/ 

DCL < C S * S A • SS f SC ) CHAR<*> VAR ♦ 

SC — • 5 s 

I START =1 NOEX< C St SA > 4-LENGTH (SA I ; 

I STOP~ INDE X{ SUBSTR (CS • 1ST ART ) *SB> ; 

IF I START>OE-I ST0P>1 THEN SC=SUBSTR ( CS • 1ST ART • ISTOP-1 ) ? 

END SPAN! 


00000100 

00000200 

00000300 

00000400 

00000500 

00000600 

00000700 

00000800 

00000900 

OOOOIOOO 

00001100 

00001200 

00001300 

00001400 


* process; 

TR IM : PROC < STR I NG) ; 

/* ALL TRAILING BLANKS ARE DELETED FROM THE STRING CS */ 

DCL STRING CHAR<*> VAR; 

DO l=LENGTH( STRING! TO 1 BY -1 WHILE ( SUBSTR ( STRING* I • 1 1 = • 

end; 

if 1=0 then string=««; else string=substr(string*i*ii; 

end trim; 


oooootoo 

00000200 

00000300 

00000400 

• 1 ;ooooosoo 
00000600 
00000700 
00000800 
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APPENDIX B 



SUBROUTINE TAYLOR < TI ,Yl .TF .YF.NCF.NEQ, 1 NIT , EPS, RANGE ) 

000001 00 


IMPLICIT REALMS (A-H.O-Z) 

00000200 


COMMON /COUNT/ KEV.KIS.IER 

00000300 


COMMON /PAR/ HMIN 

00000400 


REAL* 8 YI ( NCF , 1 > , YF<NCF, 1 ) 

00000500 


LOGIC AL*4 INI T 

00000600 


IF <.NOT.INIT> GO TO 11 

00000700 

C**** 

INITIALIZATION 

00000800 


IER = 0 

00000900 


KEY = 0 

00001000 


KIS = 0 

00001100 


NCL = NCF - 1 

00001200 


EX = 1 * D 0 / NCL 

00001300 


NEQ1 - NEQ + 1 

00001400 


CR » OSIGN(l *00* RANGE ) 

00001500 


CALL INITAL<TI *YF*NCF| 

00001600 


CALL COEFF <YF,NCF> 

00001700 


DO 18 K * 1 » NEQ 

00001800 


DO 18 J = 1 t NCF 

00001900 

18 

Y I ( J|K) = YF(J,K) 

00002000 

c**** 

COMPUTE R < H > 

00002100 

1 1 

R as 0.00 

00002200 


DO 10 K = l.NEQ 

00002300 


RJ » YKNCFtKl 

00002400 


IF (YI ( 1 « K ) . NE . 0 • DO ) RJ = PJ/YI(1.K> 

00002500 

10 

R ss DMAX1 ( R.DABSI RJ) ) 

00002600 

c**** 

COMPUTE H 

00002700 


I F (R.NE.O.DO) GO TO 23 

00002800 


H = RANGE 

00002900 


GO TO 24 

00003000 

23 

H » CR*(EPS/R)**EX 

00003100 

24 

IF {DABS(H) .GT.HMIN) GO TO 17 

00003200 

c **** 

error: h too small 

00003300 


IER = - 2 

00003400 


RETURN 

00003500 


B-l 




take a step 

00003600 

17 

DO 12 K = 1 » NEO 

00003700 

1 2 

YF(l.K) = Yl ( NCL * K 1 

00003800 


DO 13 JJ - 2.NCL 

00003900 


J s NCF - JJ 

00004000 


DO 13 K 1 ,N EQ 

000041 00 

13 

YF(l.K) = YI(J.K) ♦ H*YFI1*K> 

00004200 


TF = TI + H 

00004300 


YF( 1 #NEO 1 I - TF 

00004400 


CALL COEFF(YF,NCF) 

00004500 


KEV = KEV + 1 

00004600 


RANGE = RANGE - H 

00004700 


KIS = KIS ♦ 1 

00004800 


RETURN 

00004900 


END 

00005000 


SUBROUTINE INTERP I TI , YI * NCF , NEQ.T W 
IMPLICIT REAL *8 (A-H.O-Z) 

REAL *8 YMNCF,1)»W(1) 

H — TW - T I 
NCL = NCF - I 
DO 21 K = 1 * NEQ 

21 W(K) = YI(NCL.K) 

DO 22 JJ = 2 • NCL 
J = NCF - JJ 

DO 22 K = 1 t NEQ 

22 #(K) = YI ( J * K ) 4- H*W<K> 

RETURN 

END 


oooosioo 

00005200 

00005300 

00005400 

00005500 

00005600 

00005700 

00005800 

00005900 

00006000 

00006100 

00006200 

00006300 
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SUBROUTINE ZERO! T *F « A *TZ *B *C * N I 
IMPLICIT REAL*8< A-H,0-Z) 
DIMENSION A ( N ) *B < N ) t C ( N*N) 

C ( 1 , 1 ) =A < 1 » 

YP=-F 

B( 1) = 1 • DO/C < 1 *1) 

TZ=T+B< 1)*YP 
DO 500 K=2*N 
C(K,1 >=A< K> 

00 300 J=2*K 
C{ Ki J)=0.D0 

1 MAX— K— J-f 1 

DO 300 1 = 1 * I MAX 

300 C(K* J)=Cf K-I*J-1)*A(I > 

B( K)=O.ODO 
I MAX=K—1 

00 400 1 = 1 .IMAX 
400 B { K ) =B <K)+C{K*I)*B(I ) 
B<K)=-B(K>/C( K,K) 

YP=YP*(-F> 

TZ=TZ+B( K ) *YP 
500 CONTINUE 
RETURN 
END 


00000100 

00000200 

00000300 

00000400 

00000500 

00000600 

00000700 

00000800 

00000900 

00001000 

00001100 

00001200 

00001300 

00001400 

00001500 

00001600 

00001700 

00001800 

00001900 

00002000 

00002100 

00002200 

00002300 

00002400 
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APPENDIX C 


In the following compilation of the recurrence coeffi- 
cients for commonly used non-rational functions , it is 

I.L. 

assumed that the k— Taylor coefficients for the function 
A(t) are known and the k— coefficient (k>l) for B=f(A) is 
sought. 

A (t) = 2 a. (t-t Q ) j 
j=o J 

B (t) =2 b. (t-t ) ^ 
j=o J 

For each function, the functional relationship is listed 
first, followed by the defining differential equation and 
finally by the Taylor coefficients. Many functions such as 
sin, cos, tan are derived from a coupled differential 
system and consequently are listed together. The symbol ’ 
denotes differentiation with respect to the independent 
variable t. 


C-l 



B = exp (A) 


B’ = BA' 

b k = !V a A-i )/k 

j-1 


B = log e < A ) / a o >0 
B’ = A ’/A 

b l = V a o 

b k = (a k - Z , ja k-i b i / k) / a o . k>2 

j-x - - ” 


lo 9io (A) 

, a o >0 

aA* /A , 

a = log(e) 

aa ]/ a o 

k -1 

( “ a k • 

• = 1 ja k-j b j /k)/a o 


k>2 


B = A a _ a real, a >0 
r o 

B' = aBA’/A 

k-1 

b, = Z (a - j(a + l)/k)b.a, .)/a ft 
K j =0 3 K 3 o 


B = sin (A) , C = cos (A) , D = tan (A) , c o 4 0 
B* = CA 1 , C' - -BA 1 , D = B/C 


b, = (I ja.c, _.)/k 
* j=1 i k 3 

k 

C k = -(£ ja.b. _.)/k 
K j=l 3 k 3 


% = < b k - L^»/, 

D l 


02 



B = sinh (A) , C = cosh (A) , D = tanh (A) , c^+O 
B' - CA' , C' - BA* , D - B/C 


b k = ( j =1 ja j C k-j )/k 
k 

C k = (E j a -i b k _^)A 

* j=1 D k 3 


= (b k - 


k 

Z c . d, ) /c 
• i 3 k-j" o 
j=l J J 


The Taylor coefficients for the operations +, - , *, / are 
also included. 

C = A + B 

C k = a k + b k 


C = A - B 


c 


k 


a,_ - 


b 


k 


C = A*B 


k 

- Z 


j=o 


a, .b . 
k-J 3 


C = A/B 


c 


k 



k 

Z b. c . .)/b 

j =1 o 
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APPENDIX D 


RMAT ENTRY 

LEVEL K TYPE (C , 0 , V | E) 

# Y { 1 , 2 ) =1 . 0 $ 

12 0 Y (1, 2) 

13 5- 

1 4 -1 #2 

1 5 9 $ 

#Y (1 ,1) =Y (1) **2+3 .DO *T** 2$ 

12 0 Y (1 , 1) 

13 5 = 

14 0 Y ( 1) 

1 5 10 ** 

1 6 -1 #3 

17 1 + 

1 8 -I #4 

1 9 3 * 

1 10 0 T 

1 11 10 ** 

1 12 -1 #3 

O 1 13 9 $ 


R OP 


A(l) 


A ( 2 ) 


A (3 ) 




NASA-GSFC 


CONSTANT TABLE 


# 1=0. 5D0 # 2=1. 0D0 # 3=2 # 4=3. DO 

DMAT ENTRY 1, 2, 3, 4, 5, 6 
D MATRIX 


1, 

1 

-1 

1, 

2 

-1 

1, 

3 

-1 

1, 

4 

-1 

1, 

5 

-1 

1, 

6 

-1 

2, 

1 

-1 

2, 

2 

-1 

2 / 

3 

-1 

2 , 

4 

-1 

2, 

5 

-1 

2, 

6 

0 

3 r 

1 

-1 

3, 

2 

-ft- 

3, 

3 

-1 

3, 

4 

-1 

3, 

5 

-1 

3/ 

6 

-1 

4, 

1 

-1 

4 / 

2 

-1 

4 r 

3 

1 

4, 

4 

-1 

4, 

5 

-1 

4, 

6 

-1 

5, 

1 

-1 

5 , 

2 

1 

5, 

3 

1 

5, 

4 

1 

5, 

5 

-1 

5 , 

6 

0 

6 , 

1 

-1 

6, 

2 

1 

6, 

3 

1 

6, 

4 

1 

6 , 

5 

1 

6 / 

6 

0 














CORRESPONDENCE BETWEEN RECURRENCE MATRIX ROWS AND THE Y ARRAY 
7 1 2 23 34 45 56 61 



